Introduction
It has been well established observationally that oceanic heat flux and circulation play a significant role in determining the ice-edge location in the Bering Sea. The amount of ice transported towards the edge also affects the ice-edge location, as do wind driven and tidally induced currents. A useful way to investigate the relative role of different oceanic processes on the Bering Sea ice edge and ice drift is to carry out coupled ice–ocean simulations over several seasonal cycles. Although the simulation of all of the major effects would require a high-resolution free-surface ocean model, it is possible to gain considerable insight into the role of oceanic processes on the longer time-scale variations of the Bering Sea ice cover by constructing a coupled rigid lid ice-ocean model. This type of model also helps to identify physical processes and data needed for the simulation of ocean circulation in this region.
To achieve these goals, we have developed here a 40 km-resolution ice–ocean model of the Bering Sea. The model is very similar to the ice–ocean model of Reference Hibler,III and BryanHibler and Bryan (1987), which was applied to the Arctic, Greenland and Norwegian seas and which yielded a greatly improved ice edge in the East Greenland region. The basic difference is that the current model for the Bering Sea has a considerably higher resolution. It is also being applied to a region dominated largely by a shallow shelf in the east and a deep basin (Aleutian Basin) in the west.
By carrying out a series of sensitivity studies with this ice–ocean model, it is possible for us to examine the overall behavior of both ocean and ice in the Bering Sea. Particular emphasis is placed on determining the physical processes responsible for the seasonal and interannual ice-edge variations. The ice and ocean circulations are also examined in order to identify improvements needed in future studies.
Model Description
Basic features
The coupled ice–ocean model is described in detail in Reference Hibler,III and BryanHibler and Bryan (1987). We outline some of its basic features here. The model is constructed by coupling the two-level dynamic-thermodynamic sea-ice model of Reference Hibler,III.Hibler (1979) and Reference Hibler,III and BryanHibler and Walsh (1980) with a multi-level baroclinic ocean model (Reference BryanBryan, 1969). The sea-ice model supplies heat flux, salt flux and momentum exchange boundary conditions for the top of the ocean,whereas the ocean model gives the sea-ice model current and heat exchange information. The main forcing fields are atmospheric wind, surface air temperature and humidity. Thermodynamic heat exchanges are conserved by ensuring that any warming by the ocean in the presence of ice results in ice melting. Atmospheric heat losses and gains are calculated by the sea-ice model (even if no ice is present). In addition to thermodynamic exchange with the ice, the ocean also affects the ice drift via modification of surface currents and the oceanic tilt term in the momentum balance. The model is a robust diagnostic one with observed annual ocean temperature and salinity from Reference LevitusLevitus (1982) used to constrain the model weakly (three-year relaxation time) at all levels below the top level of the ocean, which is taken as a de facto mixed layer.
Fig. 1. Grid configuration of ice–ocean model used in the Bering Sea. Stronger relaxation (30 to 220 d) cells and open boundary cells are marked in the grid. In unmarked cells, a weak relaxation (three years) was made. The lettered grid cells denote locations of river runoff by the Yukon and Kuskokwim rivers. The grid is divided into three sectors for later time series analysis.
Fig. 2. Contours of ocean bottom topography used by the model in the Bering Sea.
Grid configuration
The grid configuration of the model, shown in Figure 1, covers the Bering Strait and encloses the Bering Sea along the Aleutian island chain. The ice portion of the model has a 40 km by 40 km horizontal resolution with a grid size of 41 by 54 in rectangular co-ordinates. The ocean portion of the model uses a spherical geometry and has a spacing of 0.36– by 0.36– in a spherical co-ordinate system whose equator is aligned with 35– Ε and passes through the North Pole. The salt flux errors caused by the geometric mismatch in the grids have been found negligible for purposes of a diagnostic simulation (Reference Hibler,III and BryanHibler and Bryan, 1987). The vertical levels are also the same as used by Reference Hibler,III and BryanHibler and Bryan (1987). Bottom topography is resolved by using differing numbers of levels at different locations. Figure 2 shows the model bottom topography contours in the analysis region. In the very shallow eastern Bering Sea only three vertical levels are considered. This treatment overestimates the depth of some of the shallower portions of this region. The grid is divided into three sectors (Fig. 1) to facilitate the interannual variability analysis.
Boundary conditions and forcing fields
The southern boundary of the model is effectively formed by the Aleutian chain, where a rigid wall is imposed for both the ice and ocean components of the model. An open boundary for the ice portion is set outside the Bering Strait (Fig. 1) whereas a wall is imposed there for the ocean portion.
A weak diagnostic relaxation to mean observed annual temperature and salinity is made with a three-year relaxation time constant in most parts of the grid, except for the two hatched regions indicated in Figure 1. In the hatched region, adjacent to the wall outside the Bering Strait, stronger relaxation with a damping constant of 30 d is imposed so that the observed temperature and salinity are closely followed near that boundary. In the hatched region, near the boundary of the western Bering Sea, stronger relaxation is also imposed with the damping constant varying from 30 d along the leftmost column to 220 d along the rightmost column of the stronger relaxation cells (Fig. 1).
For river runoff into the Bering Sea, the main source is from the Yukon River, with only minor contributions from other sources, most notably the Kuskokwim River (Reference Overland, Hood and CalderOverland, 1981). The monthly runoff rate for these two rivers from Reference Overland, Hood and CalderOverland (1981) is used. The grid cells in Figure 1 that contain letters denote cells where river runoff is specified. For the Yukon River, two cells are used to introduce its inflow.
In the standard model run no precipitation is assumed. However, in a sensitivity study, a precipitation value of 1.0 ma−1 is used throughout the whole region to estimate roughly the role precipitation plays in ice extent. This value is an approximate average of the observed data reported by Reference Overland, Hood and CalderOverland (1981).
Atmospheric forcing fields to drive the ice–ocean model are identical to those used by Reference Hibler,III, Ranelli, Johannessen, Haugan and OwensHibler and Ranelli (1990). Specifically, surface atmospheric pressures and air temperatures were interpolated from 5° by 5° latitude-longitude fields provided by John Walsh (personal communication). Daily geostrophic winds were computed from the pressure data, which were obtained by merging
Arctic data-buoy pressure analysis with NCAR daily surface analyses. The air temperatures were from monthly mean temperature fields obtained from NASA analyses. A constant relative humidity of 0.9 was assumed. The shortwave and long-wave radiations were computed using the formulation of Reference Parkinson and WashingtonParkinson and Washington (1979).
Simulation Results and Discussion
To obtain interannual results, the ice-ocean model is integrated for a three-year period (1981–83) using interannual forcing fields after an initial two-year spin-up with cyclically repeated forcing from 1980. A distorted time-stepping procedure is used with an integration time interval of 6 h for the ocean conservation equations and the ice model, and a time interval of 18 min for the ocean dynamic equations. In the model momentum equations, a constant horizontal eddy viscosity of Am = 2 × 108cm2s−1 and a constant vertical viscosity of K m = 1.0cm2s−1 are used. The horizontal value is a factor of five smaller than that used by Reference Hibler,III and BryanHibler and Bryan (1987). This reduction was made possible by the higher resolution. In the conservation equations for salt and heat, a constant lateral diffusion constant of A h= 107cm2s−1 was used; a value identical to that used by Reference Hibler,III and BryanHibler and Bryan (1987).
To approximate the effect of penetrative convection on the ice margin, two interannual simulations were carried out with different vertical diffusion constants. In an initial run, a vertical kinematic diffusion constant of K h = 0.1 cm2 s−1 was selected. This value is identical to that used by Reference Hibler,III and BryanHibler and Bryan (1987), and basically yields a model where large-scale convective overturning tends to supply more heat flux than diffusion. For a second study, a value of K h = 10.0 cm2 s−1 was selected. The idea here is to approximate the effects of penetrative convection which are not explicitly treated in the model. For comparison, a third interannual simulation using an ice-only model with a motionless fixed depth boundary layer (with no oceanic heat flux) was carried out. For analysis purposes, two other short sensitivity studies were made. They consist of a two-year simulation (1980 and 1981 atmospheric forcing) with precipitation added as rain in an uniform spatial and temporal manner; and a one-year barotropic ocean model using 1981 forcing. In the barotropic model, the stress for the ice-only case was used to drive the ocean model and the temperature and salinity fields were forced to be spatially uniform.
Ice and ocean circulation characteristics
To examine the characteristics of ocean and ice velocities, we will focus on the ice–ocean model with high diffusion and without rainfall. The ice and ocean velocities of the low-diffusion simulation were found to be similar, with the main difference being the ice-edge characteristics. The average 1981 winter surface currents for the baroclinic and barotropic ice–ocean simulations are compared in parts a and b of Figure 3. As can be seen from this figure, there is a baroclinic counter-clockwise gyre in the Aleutian Basin with some current intensification along the coast of the Soviet Union. On the Bering Sea shelf the flow is more variable and largely barotropic. In summer, the current structure over the Bering shelf was found to weaken considerably in response to weaker winds. However, the surface currents of the Aleutian Gyre show less change. Typical maximum current speeds are seen to be of the order of 5cms−1. Although not shown, there is considerable interannual variability in the surface currents in winter due to variations in the wind stress.
Fig. 3. Simulated 1981 winter (January, February and March) average ocean surface currents (from the 2nd level in the model) for: (a) full ice–ocean model and (b) barotropic ice–ocean model.
Comparison of the simulated currents with observations (Reference Kinder, Schumacher, Hood and CalderKinder, 1981; Reference Kinder, Hood and CalderKinder and Schumacher, 1981) indicates that the Aleutian Gyre circulation, and the northwestern shelf flow along the coast of Alaska, are qualitatively reasonable. However, the simulated south-eastern flow along the shelf break is unrealistic. There are several likely reasons for this. One possibility is not having a net transport of approximately 1 Sverdrup through the Bering Strait into the Arctic. Since the simulated, vertically integrated flow over the shelf is about 1 Sverdrup, the Bering Strait would supply an outlet for the northwestern flow along the coast without it having to return southeast along the shelf break. Another difficulty is the lack of strong density gradients along the shelf break in the initial observed data to which the model is diagnostically forced. Such density gradients could produce a northwestern surface current along the shelf break even though the vertically intergrated flow is in the other direction. Finally, there is an inadequate resolution of the shelf topography in this model (it is essentially treated as a flat bottom) due to limited vertical resolution.
Fig. 4. Simulated winter average ice velocity distribution for: (a) 1981 and (b) 1982.
The average winter ice velocities shown in Figure 4 are substantially larger in magnitude and simpler in pattern than ocean surface currents. Overall, the ice drift tends to be much more uniform in direction. The 1981 winter average shows a consistent southwestward drift. However, in response to wind changes, the 1982 winter ice drift tends to be more westward, at least south of St. Lawrence Island. The 1983 winter ice drift was found to be very similar to 1981. Note that there is little similarity of ice drift to the ocean currents due to the ice drift being about five times larger in magnitude than the surface currents.
That the ice drift is mainly wind driven and reasonable in magnitude is verified in Figure 5. This shows a comparison between simulated and observed ice drift over a short period from 12 to 24 February 1983. This figure shows an observed ice-buoy trackline and simulated ice-drift trackline obtained by integration of daily ice velocities from the ice-only model and ice-ocean model, respectively. The observed trackline is from Michael and Pease (1985). Note that the inclusion of ocean currents improves the agreement with the observed ice drift, although the effect of the currents is not particularly large.
Fig. 5. Comparison of observed ice-buoy drift and simulated ice drift by the ice-only model and the higher diffusion ice-ocean model from 12 to 24 February in 1983.
Ice-edge variations
Of particular interest to this study is the effect of ocean circulation on the seasonal and interannual ice-edge fluctuations in the Bering Sea. The overall seasonal behavior of the two ice–ocean models and the ice-only model are shown in Figure 6, where simulated 0.2 concentration limits are compared with observed ice limits for three typical months in 1981. Also shown in part (d) of this figure is the effect of precipitation on the ice margin. As can be seen, the ice-only model produces excessive ice, particularly in the western Bering Sea, throughout the whole time period. This excessive ice production in the west is improved substantially by the high-diffusion ice-ocean model. However, this improvement is partially at the expense of the ice extent in the eastern Bering Sea, where the maximum extent is underestimated by the high diffusion model. Also, in the eastern Bering Sea region, the rate of advance of the simulated ice edge is too slow and the retreat too fast.
Some of these deficiencies and characteristics of the seasonal cycle of ice extent are shown more clearly in Figure 7, where 3-year average seasonal ice coverages in the three sectors are illustrated. In sector I (eastern Bering Sea) we see that all models, with the exception of the ice-only model, underestimated the average seasonal cycle of ice extent. In the other sectors the reverse situation applied for all but the high-diffusion model, with excessive ice being simulated. Moreover, in sector III, even the high diffusion model yields somewhat excessive ice. Note that the best fit between observed and simulated seasonal ice
Fig. 6. Observed and simulated monthly average ice-edge locations in 1981 for (a) February; (b) May; (c) December. Part (d) of this figure shows the ice edges for the low diffusion model with surface precipitation included for February 1981.
Fig. 7. Three-year average ice coverage in three sectors in the Bering Sea. The unit is fraction of coverage of the whole sector.
extent occurs in sector II, when the high diffusion ice–ocean model is used.
The magnitude of the additional heat flux produced by the higher vertical diffusion is illustrated in Figure 8. This shows average heat fluxes in winter 1981 for the two ice-ocean models. This heat flux is the heat transfer from the deeper ocean to the top de facto mixed layer of the ocean. As can be seen, the vertical diffusion is particularly effective in the western Bering Sea near the Soviet coast, where the heat flux magnitude is as large as 120Wm−2. This value is more than doubled over the low diffusion results. There is also a greater heat flux over the shelf in the high diffusion case, although in both models the heat fluxes over the shelf are substantially smaller than over the deep water. Our conclusion is that, in these high heat flux regions, a weak stratification prevents convective overturning from bringing up the substantial amount of deeper warm water available. However, some type of penetrative convection, crudely parameterized here by a high vertical diffusion, should enable the heat to be made available to the ice.
Whereas sectors II and III have too much ice, the opposite situation applies in sector I, where too weak a
stratification allows too much convective overturning to occur. Note from Figure 8 that in this region there is less difference between the heat flux from the two ice-ocean models. As shown by Figure 6(d), the situation in this sector can be improved substantially by the inclusion of precipitation, which is known to be particularly heavy in the southeastern part of the Bering Sea. With more precipitation, a large amount of fresh water is distributed on the top of the ocean, forming a more stratified ocean. In the case of the low-diffusion ice-ocean model, this stratification reduces the convective overturning. Although it is not shown in these figures, we note that when a high diffusion term is used the effect of the rainfall is less noticeable. This is because the higher vertical diffusion causes mixing to occur relatively independent of stratification.
Fig. 8. Simulated 1981 winter average heat flux from the deeper ocean to the top 30 m de facto mixed layer. The contour interval is 15 Wm−2. (a) High diffusion ice–ocean model; (b) Low diffusion ice–ocean model.
Fig. 9. Deviation of the monthly ice coverages from the three-year monthly mean in the three sectors from 1981 to 1983. Only the months when there was significant ice (November to April) are plotted. The months 1 to 18 correspond to January to April 1981, November 1981 to April 1982, November 1982 to April 1983 and November to December 1983. The solid line represents observed data and the dashed line simulated results. The unit is fraction of ice coverage of the whole sector.
Table 1. Correlation coefficients between simulated monthly and observed interannual ice coverages shown in Figure 9.
The interannual characteristics of the ice-edge simulations in the various sectors are illustrated in Figure 9. In this figure we plot the deviation from the average seasonal cycle of ice extent in each sector for the months November through April. A more quantitative measure of the performance of different models is given by correlation coefficients in Table 1 between the simulated
interannual ice coverages and the observed variations shown in Figure 9. As can be seen from the figure and the table, the high-diffusion ice-ocean model performs best in that it gives the highest correlation in all sectors. Moreover, in sectors II and III, the magnitude as well as the variability of the interannual variations are well simulated. This can, for example, be contrasted with the ice-only model in sector III where there is a good correlation, but much too small a magnitude of simulated variability. In sector I the high diffusion model still performs best, but the correlation is weak and the magnitude of the variability is underestimated.
To summarize the ice-edge results, the high diffusion model simulates the ice-edge interannual variability in sectors II and III very well. It also simulates the seasonal cycle in these sectors reasonably well. However, no model performs well in sector I. One problem in sector I would appear to be inadequate stratification due to neglect of the large amount of precipitation occurring in that region. The difficulties in sector I emphasize the fact that simply adding diffusion everywhere is a poor substitute for including a more realistic treatment of penetrative convection via some type of boundary layer formulation. In particular, specifying a higher diffusion causes a greater mixing, regardless of the state of stratification or the amount of buoyant or mechanical energy available for deepening the mixed layer. The inclusion of a more realistic mixed layer clearly has a high priority in future model studies.
Conclusions
The main purpose of this study has been to develop a coupled diagnostic ice–ocean model of the Bering Sea to assess the effects of oceanic circulation on ice-edge variations. The results have not been totally satisfactory. In particular, whereas the coupled ice-ocean model does give significantly better results than the ice-only model for the ice margin, the ocean currents over the shelf and shelf break are not realistic. Moreover, there is not enough ice simulated in the eastern Bering Sea. The main problem with the currents is an unrealistic simulated southeastern flow along the shelf break between the counter-clockwise baroclinic currents in the Aleutian Basin and the counterclockwise barotropic flow over the shelf. The most likely cause of this unrealistic current on the shelf break is the lack of strong density gradients in this region in the “observed” data. The winter average ice velocity distributions do, however, appear reasonable and consistent with winds. A short period of integration of ice drift is in agreement with observed ice-buoy drift.
In an attempt to parameterize the effects of penetrative convection on the ice edge, two different diffusion constants were used. These two simulations result in significant differences in ice extent due to the different amounts of oceanic heat flux into the upper layer. Overall, the high diffusion works best because it suppresses the excessive ice in the western Bering Sea. On the other hand, in the eastern Bering Sea, the best results are obtained with precipitation and a low diffusion. Without this large precipitation, and the concomitant increased stratification in the ocean, there is at least a small amount of oceanic heat flux which prevents the ice from forming as far south as is observed. Clearly this combination of results amplifies the need for the inclusion of a more complete boundary layer formulation with greater vertical resolution, in which vertical mixing is made a function of stratification and turbulent kinetic energy rather than taken as a constant. The results also demonstrate the importance of precipitation (or some other fresh-water source) for simulating the location of the sea-ice margin, especially in the eastern part of the Bering Sea.
Acknowledgements
This work was funded by Greenhorne and O'Mara, Inc., via a grant from the Mineral Management Service. Mr Peter Ranelli gave us considerable help in the numerical simulations described here.
