A coupled sea ice–ocean model is developed to quantify the sea ice response to changes in atmospheric and oceanic forcing in the Bering Sea over the period 1970–2008. The model captures much of the observed spatiotemporal variability of sea ice and sea surface temperature (SST) and the basic features of the upper-ocean circulation in the Bering Sea. Model results suggest that tides affect the spatial redistribution of ice mass by up to 0.1 m or 15% in the central-eastern Bering Sea by modifying ice motion and deformation and ocean flows. The considerable interannual variability in the pattern and strength of winter northeasterly winds leads to southwestward ice mass advection during January–May, ranging from 0.9 × 1012 m3 in 1996 to 1.8 × 1012 m3 in 1976 and averaging 1.4 × 1012 m3, which is almost twice the January–May mean total ice volume in the Bering Sea. The large-scale southward ice mass advection is constrained by warm surface waters in the south that melt 1.5 × 1012 m3 of ice in mainly the ice-edge areas during January–May, with substantial interannual variability ranging from 0.94 × 1012 m3 in 1996 to 2.0 × 1012 m3 in 1976. Ice mass advection processes also enhance thermodynamic ice growth in the northern Bering Sea by increasing areas of open water and thin ice. Ice growth during January–May is 0.90 × 1012 m3 in 1996 and 2.1 × 1012 m3 in 1976, averaging 1.3 × 1012 m3 over 1970–2008. Thus, the substantial interannual variability of the Bering Sea ice cover is dominated by changes in the wind-driven ice mass advection and the ocean thermal front at the ice edge. The observed ecological regime shifts in the Bering Sea occurred with significant changes in sea ice, surface air temperature, and SST, which in turn are correlated with the Pacific decadal oscillation over 1970–2008 but not with other climate indices: Arctic Oscillation, North Pacific index, and El Niño–Southern Oscillation. This indicates that the PDO index may most effectively explain the regime shifts in the Bering Sea.
The Bering Sea supports one of the world’s most productive and valuable fisheries, immense populations of marine birds and mammals, and the subsistence activities of native communities (ARCUS 2004). This great biological productivity takes place in a dynamic ocean with substantial seasonal ice cover. Significant changes in the Bering ecosystem have been measured in recent decades. Observations indicate that large-scale ecological “regime shifts” occurred in the Bering Sea during 1976–77, 1988–89, and 1998, affecting the mix and abundance of suites of coexisting species from primary producers to apex predators (e.g., Stockwell et al. 2001; Benson and Trites 2002; Hunt et al. 2002; Hunt and Stabeno 2002). A major ecosystem shift has also been observed in the northern Bering Sea (Grebmeier et al. 2006).
Although our current understanding of the ecological regime shifts in the Bering Sea is not deterministic (Overland et al. 2008), it appears the shifts are closely linked to changes in the Bering climate. Since 1976, surface air temperature (SAT) has generally been increasing, particularly in summer; since 1996, spring consistently arrives earlier. Since 2000, warm SAT anomalies generally extend from February through November (Overland and Stabeno 2004). Long-term moored oceanographic data from the eastern Bering Sea shelf indicate that the depth-averaged spring and summer water temperatures are ∼2°C warmer in the 2001–03 mean than in the 1995–97 mean (Stabeno et al. 2002; Overland and Stabeno 2004). This warming continues until at least 2005 (Stabeno et al. 2007). However, this is not an uninterrupted warming trend; data show that the Bering Sea was very cold in 2008 (see online at http://bsierp.nprb.org/results/headlines.html).
Changes in the Bering climate in the past decades appear to be part of a larger North Pacific shift associated with warmer water temperatures and lower sea level pressures (SLPs) in the winter Aleutian low. Various modes of climate variability have been used to describe changes in Bering climate such as the Pacific decadal oscillation (PDO; Mantua et al. 1997), the North Pacific index (NPI; Trenberth and Hurrell 1994), and the Arctic Oscillation (AO; Thompson and Wallace 1998; Overland et al. 1999). The NPI provides a measure of the intensity of the mean wintertime Aleutian low, whereas the AO is the leading mode of arctic SLP variability that is also linked to the Aleutian low (Stabeno and Overland 2001). The PDO, the leading mode of sea surface temperature (SST) variability in the North Pacific, was in the positive phase from 1977 to the late 1980s, with warmer sea temperatures along the west coast of North America. Recently, the PDO has shifted between negative and positive phases without a persistent pattern. The El Niño–Southern Oscillation (ENSO) events may also influence the eastern Bering Sea through teleconnections (Niebauer and Day 1989; Overland et al. 2001).
A key feature of the Bering marine environment is the sea ice that varies substantially on multiple temporal scales in response to climate change (Niebauer et al. 1999). Winter southwestward ice advection driven by northerly/northeasterly winds and ice melt (IM) at the thermodynamic limit (warm waters in the south) are the principal processes influencing sea ice production/extent in the Bering Sea (Pease 1980). The eastern Bering ecosystem responds strongly to changes in the location of the ice edge and ocean fronts (for an overview, see Macklin et al. 2002). Ice cover affects the availability of light to the water column and the air–sea exchange of biogenic gases. Ice provides both a buoyancy source and sink at the ocean surface, and influences upwelling processes and wind-driven ocean mixing. Sea ice appears also to be important in the formation and maintenance of “cold pools” (Stabeno et al. 2002), regions of cold water in contact with the seafloor that are often associated with high nutrient distributions (Whitledge et al. 1988). The productive “Green Belt” along the Bering Sea slope region, which reflects slope upwelling processes and the shelf–basin exchange of physical and biological properties through currents, eddies, and meanders (Springer et al. 1996; Okkonen et al. 2004), may also be modified by the ice cover. The ice itself is a habitat for a vast range of organisms from ice algae to large marine mammals. Therefore, variability and change in the extent, thickness, and duration of the Bering Sea ice cover will have a major impact on ecosystem functioning and biogeochemical cycling (e.g., Grebmeier and McRoy 1989; Hunt et al. 2002). Thus, it is essential to further our understanding of the behavior of the Bering Sea ice cover in response to climate variability, through field monitoring and numerical model simulations.
A variety of numerical models has been used successfully to study sea ice in the Bering Sea. Some of the models are stand-alone sea ice models (e.g., Pease and Overland 1989). Kantha and Mellor (1989) were the first to use a coupled ice–ocean model to study the Bering Sea marginal ice zone, followed by Pritchard et al. (1990), who aimed to forecast the Bering Sea ice edge. Zhang and Hibler (1991) studied the effect of ocean circulation on seasonal and interannual ice-edge variations in the Bering Sea. Recently, Wang et al. (2009) examined seasonal variations of the Bering Sea ice cover with seasonally varying atmospheric and oceanic conditions. A sea ice model has been coupled to the Regional Ocean Modeling System (P. Budgell 2007, personal communication; available online at http://www.myroms.org/), and the coupled ice–ocean model has been used to simulate sea ice in the Bering Sea (A. Hermann 2008, personal communication; available online at http://www.ecofoci.noaa.gov/efoci_icemodels.shtml). Although these models yield considerable insights about the behavior of sea ice in the Bering Sea, they have not yet been used to quantify the sea ice response to changes in atmospheric and oceanic forcing on interannual and decadal time scales.
We have developed a coupled ice–ocean model, which we refer to as the Bering Ecosystem Study Ice–Ocean Modeling and Assimilation System (BESTMAS). In this paper, we use this model to quantify the interannual and decadal variability of the Bering Sea ice cover over the period 1970–2008. Specifically, we seek to quantify the following: 1) the mechanisms controlling the observed variability of sea ice on interannual and decadal time scales, 2) the effects of the interactions between dynamic and thermodynamic sea ice processes, and 3) the sea ice response to changes in the atmospheric and oceanic forcing. The BESTMAS model is described in section 2. In section 3, we verify the model simulation by comparison with available observations. In section 4, we quantify the various causes of the interannual variability of the Bering Sea ice cover. Concluding remarks are given in section 5.
2. Model description
a. Model components
BESTMAS is based on the coupled parallel ocean and sea ice model (POIM) of Zhang and Rothrock (2003). The sea ice model is the multicategory thickness and enthalpy distribution (TED) sea ice model (Zhang and Rothrock 2001; Hibler 1980). It employs a teardrop viscous-plastic rheology (Zhang and Rothrock 2005), a mechanical redistribution function for ice ridging (Thorndike et al. 1975; Hibler 1980), and a line successive relaxation (LSR) dynamics model to solve the ice momentum equation (Zhang and Hibler 1997). The heat equation is solved over each ice thickness category using Winton’s (2000) three-layer thermodynamic model. The TED ice model also includes a snow thickness distribution model following Flato and Hibler (1995).
The ocean model is based on the Parallel Ocean Program (POP) developed at Los Alamos National Laboratory. The POP model is an advanced version of the Bryan–Cox–Semtner-type ocean model (Bryan 1969; Cox 1984; Semtner 1986). POP includes many improvements such as incorporation of the K-profile parameterization (KPP) scheme for vertical mixing (Large et al. 1994), implementation of an implicit free-surface formulation of the barotropic mode, and adaptation to shared-memory or distributed-memory parallel computers (Smith et al. 1992; Dukowicz and Smith 1994).
Given that tidal energy accounts for 60%–90% of the total horizontal kinetic energy over the southeastern shelf region of the Bering Sea (Kinder and Schumacher 1981), we have incorporated into the POP ocean model the tidal forcing arising from the eight primary constituents (M2, S2, N2, K2, K1, O1, P1, and Q1; Gill 1982). The tidal forcing consists of a tide generating potential with corrections due to both the earth tide and self-attraction and loading following Marchuk and Kagan (1989).
b. Model grid
The model domain of BESTMAS covers the Northern Hemisphere north of 39°N. The BESTMAS finite-difference grid is based on a generalized orthogonal curvilinear coordinate system with a horizontal dimension of 600 × 300 grid points (Fig. 1). The “north pole” of the model grid is placed in Alaska. Thus, BESTMAS has its highest horizontal resolution along the Alaskan coast and in the eastern Bering Sea, with grid spacing ranging from 2 km along the Alaskan coast to 12 km along the Aleutian Chain, with an average of about 7 km for the whole Bering Sea, defined here as the area north of the Aleutian Chain up to Bering Strait (BS). There are 20 grid cells across Bering Strait, which allows a good connection between the Bering Sea and the Arctic Ocean.
The TED sea ice model has eight categories each for ice thickness, ice enthalpy, and snow depth. The centers of the eight ice thickness categories are 0, 0.38, 1.30, 3.07, 5.97, 10.24, 16.02, and 23.41 m. The POP ocean model has 30 vertical levels of varying thicknesses to resolve surface layers and bottom topography. The first 13 levels are in the upper 100 m, and the upper 6 levels are each 5 m thick. The model bathymetry (Fig. 1b) is obtained by merging the International Bathymetric Chart of the Arctic Ocean (IBCAO) and 5-minute gridded elevations/bathymetry for the world (ETOPO5) datasets (see Holland 2000).
c. Forcing and lateral boundary conditions
BESTMAS is forced by daily National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis (Kalnay et al. 1996) surface forcing fields from 1970 to 2008. The reanalysis forcing consists of surface winds, SAT, specific humidity, precipitation, evaporation, downwelling longwave radiation, and cloud fraction. SAT and cloud fraction are used to calculate downwelling shortwave radiation following Parkinson and Washington (1979). Model forcing also includes river runoff of freshwater in the Bering and Arctic Seas. For the Bering Sea, monthly climatological runoffs of the Anadyr, Yukon, and Kuskokwim Rivers are used (see Table 1 for sources). For the Arctic Ocean, monthly climatological runoffs of the Pechora, Ob, Yenisei, Olenek, Yana, Indigirka, Kolyma, Mackenzie, Dvina, Lena, Khatanga, Taimyra, and Piasina Rivers are from the Alfred Wegener Institute (Prange and Lohmann 2004) based on data from the Global River Data Center in Koblenz, Germany (see also online at http://efdl.cims.nyu.edu/project_aomip/forcing_data/rivers/overview_awi.html).
Although BESTMAS has a large model domain that includes the Arctic and the North Pacific, realistic lateral open boundary conditions are still necessary to create the right water masses and transports in the Bering Sea. We have modified the POP ocean model to incorporate open boundary conditions so that BESTMAS is able to be one-way nested to a lower-resolution but global POIM (Zhang 2005). Monthly-mean open boundary conditions of ocean temperature, salinity, and sea surface height from the global POIM are imposed at the southern boundaries along 39°N.
3. Model results and observations
The BESTMAS model is initialized from the global POIM sea ice and ocean conditions of 1 January 1969 (Zhang 2005). BESTMAS is then integrated from 1969 to 2008 using the daily NCEP–NCAR reanalysis forcing. Because we consider the first year as still in spinup mode, we present here results for the Bering Sea over the period 1970–2008. We evaluate model results by comparisons to available observations.
a. Sea ice extent and thickness and model–data comparisons
Figure 2 shows the simulated and observed seasonal and interannual variations of monthly-mean sea ice extent over 1970–2008. The simulated seasonal cycle of ice coverage is generally in phase with the observed cycle, with a few exceptions: model-simulated ice advance in the fall (October–December) and winter (January–March) of 1984 and 1985 is about one month later than observed. The model agrees with the observed winter ice extent maximum, with a 2% bias in March ice extent (Fig. 3a). However, it tends to underestimate the onset of ice in November–January and then overestimate ice extent during April–May (Fig. 3a), mainly in the western Bering Sea (Fig. 4c). Despite this, the interannual variability is reasonably close to observations; the monthly model–data correlation from January to May over 1970–2008 ranges from 0.6 in January to 0.8 in May (Fig. 3a).
The simulated 1970–2008 mean ice thickness averaged over the entire ice extent in the Bering Sea is generally below 1.2 m (Fig. 3b), and the simulated local ice thickness is generally less than 2.5 m (Fig. 4). The simulated sea ice is significantly thinner along the southern coasts of the St. Lawrence (63.42°N, 170.40°W) and St. Matthew (60.42°N, 172.75°W) Islands than in the surrounding areas (Fig. 4), indicating the effects of winter polynya formations that are often observed in those coastal areas (e.g., Stringer and Grove 1991; Grebmeier and Cooper 1995). There are very few ice thickness measurements in the Bering Sea available for model validation. Two moored upward-looking sonars (ULSs) were deployed south of St. Lawrence Island in 1999 as part of a winter polynya field experiment (Drucker et al. 2003), which yielded ∼4-month records of ice thickness. Figure 5 shows a comparison of the model-simulated daily-mean ice thicknesses and the ULS-derived estimates at those two sites (ULS-15 and ULS-42) in the winter of 1999. The model catches some of the ULS-observed ice thickness peaks (e.g., around days 23 and 63) but misses some of the other peaks (e.g., days 29 and 84). The model simulates two peaks between days 42 and 54 that occur 1 or 2 days earlier than the observations. Both model results and ULS measurements show considerable daily variability in ice thickness, reflecting the dynamic nature of the Bering Sea ice cover. The sharp increase in ULS-measured ice thickness on days 44 and 52 is likely a buildup of ridged ice (Drucker et al. 2003). BESTMAS simulates ice ridging and reproduces these two observed ice thickness peaks with some success, suggesting that it is important to incorporate ice ridging processes in models to capture the sometimes large change in ice thickness on daily time scales.
b. Ocean circulation and transports, tidal effects, and model–data comparisons
The general pattern of Bering Sea ocean flows is comparatively well known from observations (Stabeno et al. 1999). The BESTMAS model captures basic features of the observed major currents in the Bering Sea, such as the Alaskan Stream, the Aleutian North Slope Flow, the Bering Slope Current, the Kamchatka Current, and the on-shelf flows (Fig. 6; see also Fig. 2 in Stabeno et al. 1999). The model has inflows through the Aleutian passes and outflow through Bering Strait (described later). The model shows considerable seasonal variability in the strength of the 15-m ocean circulation (Figs. 6a,b), with all the major currents being significantly stronger in winter than in summer. This seasonal change in ocean circulation is generally attributed to a seasonally variable Pacific–Arctic SLP head and seasonally varying wind circulation (Wang et al. 2009).
The incorporation of tidal forcing (section 2a) in the model is found to modify the simulation of ocean flows through tidal rectification. For example, tidal forcing induces a strong clockwise circulation around the Pribilof Islands (Fig. 7a), which is qualitatively in agreement with Kowalik and Stabeno (1999) and Hermann et al. (2002). Without tidal forcing, such a clockwise circulation is not simulated (Fig. 7b). Tidal forcing also influences the simulation of ice thickness. The difference in simulated ice thickness with or without incorporating tidal forcing is often up to 0.1 m or 15% over a large area in the Bering Sea (Fig. 8). The difference is particularly noticeable in the central-eastern Bering Sea, where tides are strong (Kinder and Schumacher 1981). This is because tides affect ice motion and deformation (Pearson et al. 1981), and therefore spatial redistribution of ice mass. Tides also alter the spatial redistribution of ice mass by modifying ocean flows (Fig. 7). For example, by examining the difference in the ocean surface velocity between the simulations with and without tidal forcing, we have found that the velocity difference vectors tend to drive ice into the area north of St. Matthew Island and out of the area southwest of St. Lawrence Island. This leads to thicker ice north of St. Matthew Island and thinner ice southwest of St. Lawrence Island (Fig. 8).
A transect roughly parallel to the shelf break in the central Bering Sea (transect A in Fig. 1) is used to examine the effect of on-shelf oceanic volume and heat transports on changes in sea ice in the Bering Sea. Because river water mass is not added into the model domain but is simulated by changing ocean salinity in the model, the simulated net volume transport across transect A must, by mass conservation, balance the outflow at Bering Strait (Fig. 9a). The simulated 1970–2008 annual-mean outflow at Bering Strait is 0.75 Sv (1 Sv = 106 m3 s−1), in good agreement with the observational estimate of about 0.8 Sv (Roach et al. 1995). Observations of the Bering Strait outflow show significant interannual variability (Woodgate et al. 2005, 2006), which is generally well captured by the model, at least within observational error bars (Fig. 9a).
Like the vertically integrated volume transport, the vertically integrated on-shelf heat transport across transect A (HTA; calculated relative to −1.8°C) also shows significant interannual variability (Fig. 9b). The January–May mean on-shelf heat transport is significantly correlated with the corresponding volume transport (R = 0.64). Not surprisingly, the heat transport is also significantly correlated with the simulated SST over the whole Bering Sea and with the PDO that represents a significant portion of SST variability over the whole North Pacific (Table 2). However, both the simulated heat transport at transect A and SST in the whole Bering Sea are not significantly correlated (R < 0.4) with some other climate indices such as AO, NPI, and ENSO over 1970–2008; therefore, the corresponding correlation values are not listed in Table 2.
c. SST, SAT, and model–data comparisons
Figure 10a compares the simulated May SST in the southeastern Bering Sea with that from NCEP–NCAR reanalysis. The NCEP–NCAR reanalysis SST is derived from the global daily high-resolution Reynolds SST analyses using satellite and in situ observations (Reynolds and Marsico 1993; Reynolds et al. 2007; Kalnay et al. 1996). BESTMAS captures most of the variability in the observation-derived May SST, with a correlation of R = 0.79. The simulated annual or January–May mean SST in the whole Bering Sea is strongly correlated with the corresponding NCEP–NCAR reanalysis SAT (Figs. 10b,c; Table 2), which is not surprising because the latter is a key component of the model forcing. Both SAT and SST averaged over the whole Bering Sea show that winter 1976 is the coldest over the period 1970–2008, whereas the winters of 1979 and 1996 are the warmest (Fig. 10c). Since 2005, the Bering Sea has been cooling significantly and winter 2008 is one of the coldest, whereas winter 2007 SAT is close to the 1970–2008 average.
The different winter conditions of four sample years of 1976 (cold), 1996 (warm), 2007 (average), and 2008 (cold) are illustrated in Fig. 11, which plots fields of winter surface temperature and SAT. Here, surface temperature is the temperature on the sea ice surface or the ocean surface. Thus, the surface temperature south of the ice edge (indicated in Fig. 11) is SST. SST under the ice is close to the freezing point (−1.8°C in the model). The spatial pattern of the simulated surface temperature is generally in agreement with that of the NCEP–NCAR reanalysis surface temperature (note some irregularities in the reanalysis data for winter 1976 near 61°N, 176°E, where the reanalysis surface temperature does not match the observed ice edge; Fig. 11b). The spatial pattern of either simulated or reanalysis surface temperature is in turn close to that of the reanalysis SAT, indicating the effect of air–ice–sea interactions. However, in a large area immediately south of the ice edge, both the simulated and observed SSTs are generally above zero, even though SAT is below zero, particularly in the cold winters of 1976 and 2008. This is because of the significant heat storage in the ocean mixed layer in summer and heat supply in fall and winter to the mixed layer through vertical mixing (Stabeno et al. 2009) and northward heat transport (Clement Kinney et al. 2009; see also Fig. 9b), which is greater than the surface heat loss resulting from winter atmospheric cooling.
d. Interannual variability of sea ice
Figures 12a,b illustrate the interannual variability of the Bering Sea ice cover by showing time series of ice extent and volume averaged over January–May, a period when a significant portion of the Bering Sea is covered by ice (Fig. 3). In general, the model performs well; it overestimates the January–May mean sea ice extent by only 3% on average over the period 1970–2008 (Fig. 12a). It also captures well the variability of observed ice extent (OIE) with a model–data correlation of R = 0.78 and a root-mean-square (RMS) error of 13%. The simulated total ice volume in the Bering Sea (Fig. 12b) is similarly correlated with the observed ice extent (R = 0.76; Table 2); as expected, it is even more strongly correlated with the simulated ice extent (R = 0.96). Both ice extent and ice volume reflect the substantial interannual differences of the Bering Sea ice cover (Figs. 12a,b). Being the coldest year over 1970–2008, 1976 has the highest total ice volume and its ice cover is generally thicker over a larger area than other years such as 1996 (warm), 2007 (average), and 2008 (cold; Fig. 13). Conversely, 1996, a warm year, has the second lowest ice volume, and St. Matthew Island is almost ice free in winter (Fig. 13b). The simulated ice volume (SIV) in the Bering Sea increases gradually from 1979 to 1995 and then drops steeply in 1996. It doubles from 1996 to 2000 and then decreases significantly from 2000 to 2005. From 2005 to 2008, the ice volume and extent recover, and 2008 is one of the heaviest ice years (Figs. 12a,b, 13d).
4. Quantifying the causes of sea ice interannual variability in the Bering Sea
a. Changes in sea ice with varying SAT and SST
Bering Sea sea ice variability is due to the interplay among the climate components of the atmosphere, ocean, and sea ice itself. For example, the 2000–05 warming in the southeastern Bering Sea shelf occurred in conjunction with a mild air mass during the winter and a shorter ice season (Stabeno et al. 2007). Over the period 1970–2008, the January–May mean simulated ice volume and observed ice extent are highly negatively correlated with the corresponding reanalysis SAT (Table 2); the two time series of ice volume and SAT (also plotted in Fig. 10c) are almost symmetric horizontally (Fig. 12b). Lower SAT, often in association with lower downward longwave radiative flux (not shown), causes stronger surface cooling, leading to greater ice growth (IG). The opposite is true with higher SAT. For example, the 1976 maximum ice volume/extent occurs when both the January–May mean SAT drops to the lowest point in the record (Fig. 12b) and the winter SAT is below 0°C over almost the entire Bering Sea (Fig. 11c). Conversely, winter 1996 has a lower than usual ice volume and a higher than usual SAT (Figs. 11f, 12b). The increase (decrease) in ice volume during 2000–05 (2005–08) is associated with increasing (decreasing) SAT.
The simulated ice volume and observed ice extent are also correlated with the simulated SST in the Bering Sea and PDO, which represents a significant portion of the observed SST variability in the North Pacific (Table 2). This suggests that ocean heating may have an important role in defining ice extent and affecting ice thickness. On the other hand, the ice cover has an equally important role in modifying SST such that an ocean thermal front often forms along the ice edge, as is demonstrated by the close match between the simulated SST front and the location of the observed ice edge (Fig. 11). The simulated ice volume and observed ice extent, like the simulated on-shelf heat transport at transect A and SST averaged in the whole Bering Sea, are found to correlate poorly with other climate indices such as AO, NPI, and ENSO over 1970–2008, with correlation values generally below 0.3. This suggests that these climate indices are not particularly skillful at determining changes in the Bering Sea.
b. Changes in sea ice with varying atmospheric circulation
Wind circulation affects Bering Sea ice motion and ice mass advection. Ice motion is described by ice velocity, whereas ice mass advection is described here by ice mass convergence [−∇ · (uh)], where u is ice velocity and h is ice thickness. Thus, changes (loss or gain) in ice thickness due to ice mass advection (or ice mass convergence) quantitatively describe the effects of wind-driven ice motion and deformation on the spatial redistribution of ice mass. In winter, considerable ice loss due to ice mass advection (ILDA) generally occurs in the north and east of the Bering Sea, whereas ice gain occurs in the south and west (Fig. 14), indicating a large-scale ice movement from north and east to south and west, as described by Pease (1980). Local changes (loss or gain) in ice thickness due to ice mass advection can be up to ±0.5 m month−1. This general pattern of ice mass advection tends to expand ice extent southward significantly (Fig. 14).
Such an ice mass advection pattern is attributed to strong southwestward ice motion (Fig. 14) driven primarily by the strong northeasterly winds prevailing in winter (Pease 1980). Surface ocean circulation has a secondary role in ice mass advection because it is often in the opposite direction of ice motion (Figs. 6a, 14). The winter atmospheric circulation in the Bering Sea is often characterized by relatively high SLP in the west and north and low in the east and south because of the formation of a high pressure cell over a large area in Russia and a low pressure cell over the Aleutian Islands, leading to generally strong northeasterly winds in winter as these cells intensify (Figs. 15a–d; Pease 1980; Trenberth and Hurrell 1994; Niebauer et al. 1999; Overland et al. 1999). The strong winds, in conjunction with generally thinner ice of less strength, cause the simulated speed of winter ice motion in the Bering Sea (Fig. 14) to be greater than that in most of the Arctic Ocean (Zhang and Rothrock 2003).
One effect of the wind-driven ice mass advection is that ice is moved away from the areas downwind of the St. Lawrence and St. Matthew Islands and piled up in the areas upwind of these islands (Fig. 14; see also Pease 1987). This is illustrated by ice loss (blue color) to the south of these islands and ice gain (red color) to the north (Fig. 14). Also, ice mass advection as well as ice ridging are likely to cause the sometimes substantial daily variability in ice thickness simulated by the model and measured by the ULSs deployed in the polynya area south of St. Lawrence Island in 1999 (Fig. 5; Drucker et al. 2003).
On average over 1970–2008, the total ice loss (gain) due to ice mass advection during January–May in mainly the north and east (south and west) of the Bering Sea is about 1.4 × 1012 m3 over the 5-month period (Fig. 12c). This amount is almost twice that of the January–May mean total ice volume (0.76 × 1012 m3; Fig. 12b). In other words, on average, almost twice as much ice as the total ice volume moves from north and east to south and west of the Bering Sea. Note that the simulated 1948–2000 annual-mean sea ice export from the Arctic Ocean to the subarctic seas (mainly through Fram Strait) is about 2.7 × 1012 m3 yr−1 (Rothrock and Zhang 2005), and BESTMAS estimates that the January–May mean ice export at Fram Strait over 1970–2008 is 1.3 × 1012 m3 over the 5-month period. This means that the north/east to south/west ice mass advection during January–May in the Bering Sea is comparable in magnitude to the arctic ice export at Fram Strait.
The amount of ice moving southwestward fluctuates considerably from year to year because of the substantial interannual variability of the winter northeasterly winds, which are linked to the location and intensity of the Russian high pressure and Aleutian low pressure centers (Fig. 15; see also Overland et al. 1999). The extent and thickness of the ice cover in the southern Bering Sea are controlled by the magnitude and pattern of ice mass advection, which contributes to either an increase or a decrease of more than half of the long-term January–May mean ice volume (0.76 × 1012 m3). During the heavy ice years of 1976 and 2008 (Figs. 13a,d), ice mass advection is stronger and more expansive than usual and a large area south of 60°N gains ice due to ice mass advection (Figs. 14a,d), with a total ice gain up to 1.8 × 1012 m3 in 1976 (Fig. 12c). This is because the northeasterly winds in winter 1976 and 2008 are stronger than usual such that the wind anomaly vectors have a significant southward component over a large area, particularly south of 60°N (Figs. 15e,h). In 1996, one of the lightest ice years (Fig. 13b), the amount of ice being advected (0.98 × 1012 m3; Fig. 12c) is among the lowest and there is little ice gain in the area south of 62°N (Fig. 14b). This is because the northerly winds are so weak that the wind anomaly vectors are strongly northward (Fig. 15f). There is more ice in 2007 than in 1996, and ice mass advection in 2007 is somewhat stronger than in 1996 because of relatively stronger northerly winds or a weaker northward wind anomaly (Figs. 13c, 14c, 15g).
The close link between the wind-driven ice mass advection and ice extent/thickness is further shown by the significant correlation between the simulated ice volume and observed ice extent and the total ice loss due to ice mass advection (Fig. 12c; Table 2). The total ice loss due to ice mass advection is also correlated with the northerly wind velocities (NWV) averaged over the Bering Sea (R = 0.48; Table 2). All this suggests that atmospheric circulation, through ice mass advection processes, has a dominant role in shaping the Bering Sea ice cover. However, the correlation between the northerly wind velocities and the simulated ice volume or observed ice extent is not significant (Table 2). This suggests that the total ice loss or gain due to ice mass advection is a better proxy to describe winter sea ice response to wind forcing than wind velocities themselves.
c. Effects of the interactions among dynamic and thermodynamic sea ice processes
Ice production is a useful parameter to quantify the effects of the interactions among dynamic and thermodynamic sea ice processes. Here, ice production (Fig. 16) represents the net effects of all the air–ice–ocean thermodynamic processes modified by dynamic processes such as wind-driven ice mass advection (Zhang et al. 2008). Generally, winter ice production is strongly positive in the northern Bering Sea, indicating substantial thermodynamic ice growth that can be up to 0.5 m month−1 locally. On average over 1970–2008, the total ice growth during January–May is 1.4 × 1012 m3, close to the total ice loss or gain due to ice mass advection. The strong growth is attributed to surface atmospheric cooling with SAT much lower than 0°C (Fig. 11). The northern Bering Sea is also subject to significant ice loss due to wind-driven southward ice mass advection (Fig. 14). As ice moves southward, a greater area of open water/thin ice is created in the north, which enhances ice growth because growth rates over open water/thin ice are higher than over thick ice (e.g., Maykut 1982). Thus, the strong thermodynamic growth in the north is attributed to the combined effects of atmospheric cooling and ice opening up because of ice mass advection. This is why the January–May mean simulated total ice growth (Fig. 12d) is strongly correlated with both the January–May mean SAT (Fig. 10b; Table 2) and the total ice loss (gain) due to ice mass advection January–May (Fig. 12c; Table 2).
Southward ice mass advection also induces thermodynamic ice melt. Although winter ice growth is high in the north, significant ice melt occurs in the south, mainly near the ice edge where ice melt can be as high as −0.5 m month−1 with substantial interannual variability (Fig. 16). Given that SAT is generally below 0°C in winter even in a large area south of the ice edge (Fig. 11), it is not atmospheric heating (reflected in SAT) that causes the melting near the ice edge. Rather, it is the oceanic heat that melts the ice as winds drive ice southward into relatively warm surface waters (see also Pease 1980). Note that either the simulated or the reanalysis winter SST near the ice edge is generally above 0°C, even though SAT is below 0°C (Fig. 11). This is because of the significant heat storage in the ocean mixed layer in summer and heat supply in fall and winter to the mixed layer through vertical mixing and northward heat transport. The high (negative) correlation (Table 2) between the on-shelf heat transport across transect A and the simulated ice volume or observed ice extent indicates the importance of northward heat transport, which tends to suppress ice growth in the shelf region. The mixed layer heat contributes to strong melting in the ice-edge areas where there is also a large ice gain due to ice mass advection (Figs. 14, 16). In other words, the combination of ice mass advection and the ocean thermal front at the ice edge largely defines the location of the ice edge or the shape and size of ice extent in winter (see also Pease 1980).
In April and May, SAT increases; as a result, ice melting is more widespread and not mainly confined near the ice edge (not shown). Thus, in spring, atmospheric heating has a more prominent role in ice melting than oceanic heating. This is why the total ice melt January–May in the Bering Sea (Fig. 12d) is closely correlated with the simulated SST and reanalysis SAT (Table 2). It is even more highly correlated with the total ice loss due to ice mass advection (Table 2), indicating again the prominent role of ice mass advection in ice melting in the Bering Sea. On average over 1970–2008, the amount of total ice melt January–May is substantial at 1.5 × 1012 m3, with significant interannual fluctuations (Fig. 12d).
5. Concluding remarks
We have developed a coupled sea ice–ocean model (BESTMAS) to quantify the changes in sea ice processes in the Bering Sea in response to atmospheric and oceanic forcing on interannual and decadal time scales. Although model deficiencies exist at least in part because of uncertainties in model parameterization, forcing, and initialization, the model captures much of the observed seasonal and interannual variations of the Bering Sea ice cover. It also captures the basic features of the upper-ocean circulation and the large-scale spatiotemporal variability of SST in the Bering Sea. Model results suggest the necessity of incorporating ice ridging processes to capture the occasionally large daily change in ice thickness due to ice mass advection and deformation in the polynya area south of St. Lawrence Island. Model results also indicate that tidal forcing affects the spatial redistribution of ice mass by up to 0.1 m or 15% over a large area in the central-eastern Bering Sea by modifying ice motion and deformation and ocean flows.
The strong correlations among the simulated ice volume, observed ice extent, SAT, and SST over 1970–2008 indicate the sensitivity of sea ice to changes in atmospheric and oceanic forcing and the effect of air–ice–sea interactions. Interannual variability is considerable for the Bering Sea ice cover over 1970–2008, during which the region underwent periods of significant cooling and warming according to the reanalysis SAT and the simulated SST. The ice cover responds to the cooling and warming such that the January–May mean simulated ice volume (observed ice extent) can be as high as 1.3 × 1012 m3 (1.1 × 1012 m2) in 1976 or as low as 0.46 × 1012 m3 (0.54 × 1012 m2) in 1996. More recently, SAT–SST decreases during 1996–2000, increases during 2000–05, and decreases again during 2005–08; the ice cover changes accordingly, with a January–May mean simulated ice volume (observed ice extent) of 0.94 × 1012 m3 (0.72 × 1012 m2) in 2000, 0.52 × 1012 m3 (0.55 × 1012 m2) in 2005, and 0.97 × 1012 m3 (0.84 × 1012 m2) in 2008 (one of the heaviest ice years in the record).
As suggested by Pease (1980), changes in ice extent and thickness in the Bering Sea are most profoundly affected by wind-driven ice mass advection processes. The model allows us to quantify these processes. The Bering Sea winter features strong northeasterly winds that drive ice faster than usually observed in the Arctic Ocean, causing ice mass advection comparable in magnitude to the arctic ice export via Fram Strait. On average in the model, about 1.4 × 1012 m3 of ice, almost twice as much as the mean total ice volume in the Bering Sea (0.76 × 1012 m3), moves from the northeast to the southwest of the Bering Sea during January–May, causing significant southward expansion of ice extent. Because the atmospheric circulation varies substantially from year to year, the amount of ice moved fluctuates, from as much as 1.8 × 1012 m3 in the heavy ice year of 1976 to about half of that amount in the light ice year of 1996.
Ice mass advection also interacts with thermodynamic sea ice processes. It enhances thermodynamic ice growth in the north by increasing the areas of open water or thin ice. The average ice growth during January–May is 1.3 × 1012 m3, ranging from 0.90 × 1012 m3 in 1996 to 2.1 × 1012 m3 in 1976. Ice mass advection also enhances thermodynamic ice melting in the south by bringing ice into warm waters, causing considerable ice melting in the ice-edge areas. The average ice melt during January–May is 1.4 × 1012 m3, ranging from 0.94 × 1012 m3 in 1996 to 2.0 × 1012 m3 in 1976. In other words, ice mass advection modifies SST by cooling surface waters and forming a temperature front at the ice edge, whereas oceanic heat constrains the advection-induced southward ice expansion by melting ice near the ice edge. Thus, the interplay between ice mass advection and the ocean thermal front across the ice edge largely defines the shape and size of Bering ice extent (Pease 1980). The substantial amount of ice melting and interannual fluctuations are likely to significantly affect the meltwater distribution and ocean stratification in a large area in the southern Bering Sea, including the slope region (Green Belt) of great biological productivity.
The observed ecological regime shifts during 1976–77 and 1998 are reflected in significant changes in sea ice. The 1976–77 regime shift is associated with the heaviest ice year (1976), and the 1998 regime shift occurred with a significant increase in ice volume after two consecutive light ice years (1996–97). The correspondence between changes in sea ice and the 1988–89 regime shift is less obvious. However, the periods of 1976–77, 1998, and 1988–89, during which ecological regime shifts were observed (e.g., Stockwell et al. 2001; Benson and Trites 2002; Hunt et al. 2002; Hunt and Stabeno 2002), are all marked by significant changes in SAT, SST, and thus the PDO index (Figs. 9b, 10). Changes in sea ice, SAT, and SST all correlate significantly to the PDO but not to the AO, NPI, and ENSO over 1970–2008, indicating that the PDO may explain the regime shifts in the Bering Sea most effectively.
The results provoke speculation on a scenario of Bering Sea air–ice–ocean interactions in winter of a heavy ice year. Generally, a heavy ice year occurs with a relatively greater SLP difference between the western and eastern Bering Sea and stronger northeasterly winds, which drive more ice into the south and west of the Bering Sea, cooling the surface air and the upper ocean and melting ice there. Stronger northeasterly winds also bring more cold air into the Bering Sea from the north, thus further cooling the surface air and the upper ocean. The decreased SAT and the increased creation of open water/thin ice area due to ice mass advection enhance thermodynamic ice growth in a large area in the northern Bering Sea. In addition, because the increased presence of sea ice in the west tends to cool SAT and SST there, the relatively high SLP in the western Bering Sea and on the nearby land (Figs. 15a,d,e,h) may tend to increase further, thus likely strengthening the northerly winds. In a light ice year, a similar but reversed mechanism in the air–ice–ocean interactions in the Bering Sea would occur (Figs. 15b,f). Here, longwave radiation is not mentioned, but its effect is likely similar to that of SAT related to surface turbulent heat flux. Note that such air–ice–ocean interactions are not fully captured by BESTMAS because it is not coupled to an atmospheric model. However, some of the effects of the air–ice–ocean interactions are likely present in the NCEP–NCAR reanalysis atmospheric forcing used to drive BESTMAS.
This work is supported by NSF Grant ARC0611967. We thank R. Drucker for assistance with the ULS ice thickness data, R. Lindsay for constructive comments, and Dr. S. Okkonen for helpful suggestions.
Corresponding author address: Jinlun Zhang, Polar Science Center, Applied Physics Laboratory, College of Ocean and Fishery Sciences, University of Washington, 1013 NE 40th St., Seattle, WA 98105. Email: firstname.lastname@example.org