Multiple-Scale Variations of Sea Ice and Ocean Circulation in the Bering Sea Using Remote Sensing Observations and Numerical Modeling

Abstract

The Bering Sea is located between the Aleutian Low and Siberian High, with strong seasonal variations in the oceanic circulation and the sea ice coverage. Within such a large-scale system, the physical processes in the Bering Sea carry interannual variability. The special topography in the Bering Sea traps a strong jet along the Bering Slope, whose instability enriches the eddy activity in the region. A Regional Oceanic Modeling System (ROMS), coupled with a sea ice module, is employed to study multiple-scale variability in the sea ice and oceanic circulation in the Bering Sea for interannual, seasonal, and intra-seasonal eddy variations. The model domain covers the whole Bering Sea and a part of the Chukchi Sea and south of Aleutian Islands, with an averaged spatial resolution of 5 km. The external forcings are momentum, heat, and freshwater flux at the surface and adaptive nudging to reanalysis fields at the boundaries. The oceanic model starts in an equilibrium state from a multiple year cyclical climatology run, and then it is integrated from years 1990 through 2004. The 15 year simulation is analyzed and assessed against the observational data. The model accurately reproduces the seasonal and interannual variations in the sea ice coverage compared with the satellite-observed sea ice data from the National Snow and Ice Data Center (NSIDC). Sea surface temperature and eddy kinetic energy patterns from the ROMS agree with satellite remote sensing data. The transportation through the Bering Strait is also comparable with the estimate of mooring data. The mechanism for seasonal and interannual variation in the Bering Sea is connected to the Siberia-Aleutian index. Eddy variation along the Bering Slope is discussed. The model also simulates polynya generation and evolution around the St. Lawrence Island.

1. Introduction

The Bering Sea is a semi-enclosed sea bounded on the northwest by Siberia, on the west by the Kamchatka Peninsula, on the east by Alaska, and on the south by the Aleutian Islands. The Bering Strait, 85 km in width and 50 m in depth on average, is the only ocean gateway between the Pacific and Arctic Oceans. The prominent feature in the bottom bathymetry in the Bering Sea is the Bering Slope, which runs through the Bering Sea from the west to the east, with a sharp change in the water depth from 200 meters on the northeastern side to 2000 meters on the southwestern side within less than 100 km in width, see Figure 1. The oceanic circulation and other physical and biological processes are also significantly characterized by the presence of the Bering Slope. Due to the interaction between the baroclinic instability and the slope, numerous mesoscale eddies are generated along the Bering Slope [1,2,3,4]. Significant eddy activities along the slope partially contribute to an enhanced biological productivity zone called the “green belt” [5,6,7].

The Bering Sea is located at the interface between the north of the Aleutian Low and the southeast of the Siberian High. A large seasonal variation is observed in the Aleutian Low. In winter, the Aleutian Low is very strong, with a strong northeasterly wind over the Bering Sea, and in summer, the Aleutian Low becomes weak and even disappears, which results in a weak southwesterly wind blowing over the Bering Sea. The northeasterly wind in the winter carries the cold air from the polar area and causes the sea water to freeze and form sea ice. The southwesterly wind in summer brings warm air from the south to the Bering Sea, contributing to the sea ice melting [8]. This large seasonal variation in forcing is reflected in the oceanic circulation field and sea ice distribution [9,10,11]. The interannual change in the local winds causes the sizeable interannual variability of the Bering Strait fluxes in volume, heat, and freshwater, eventually resulting in the interannual variability of the sea ice [12,13]. Therefore, the seasonal variation of sea ice is the dominant phenomenon in the Bering Sea [14,15,16,17].

The variation in the Bering Sea’s circulation is also affected by the large-scale circulation in the Pacific Ocean. South of the Aleutian Islands, the Aleutian Stream flows westward on the northern branch of the subpolar gyre. Some significant parts of the Aleutian Stream leak into the Bering Sea and combine with local circulation [18]. The variation in the intensity of the Aleutian Stream, which carries the warm water, affects the Bering Sea circulation and sea ice variation [9,11,19,20].

Numerous observations of sea-ice coverage and hydrological variables have been conducted in the Bering Sea. The sea ice coverage is determined by means of different microwave radiation bands reflecting from the sea ice and open water, as observed from space-borne microwave radiometers [21]. This data can be traced back to the 1970s and provides information on the sea-ice coverage trend in the Bering Sea [22,23,24]. However, most available hydrological observations, such as the vertical profiles of temperature, salinity, and velocity, are located around seasonal ice-free areas, and there is a lack of observations under the sea ice that prohibits quantitative understanding of sea ice variations due to sea ice thickness and water column depth.

The past decades have witnessed increasing applications of numerical modeling to the study of sea ice and circulation variations in the Bering Sea, such as the studies by Clement et al. [25], Wang et al. [20], and Danielson et al. [9]. This ocean-only model is applied to simulate the oceanic circulation in the Bering Sea during the ice-free summertime [26,27], and cannot be used to investigate the most dominant phenomenon in the Bering Sea, the seasonal variation of the sea ice. Pritchard et al. [28] made the first effort to use an ice-ocean-coupled model, though their model’s resolution was coarse (~50 km). Clement et al. [25] used a ~9 km pan-Arctic coupled ice-ocean model and reproduced the ocean transports and interannual variability of the sea ice cover in the Bering Sea, but their model did not include the effects of tides. Wang et al. [20] developed a ~9 km coupled ice-ocean model (CIOM) based on the Princeton Ocean Model (POM) to investigate the ocean circulation and thermodynamic features in the Bering Sea. Danielson et al. [9] applied a sea ice model coupled with the ROMS implemented by Budgell [29] to study the dominant modes of interannual variability in the thermohaline and ice fields over the Bering Sea shelf. They also compared the model results with the satellites’ remote sensing and in-situ observations to examine the performance of the model. However, this study focused mainly on the vertical structure of tidal currents and the temperature and salinity fields. Our study concentrates more on the sea ice variability, ocean eddy kinetic energy, and mesoscale eddy transport.

Overall, multiple-scale variations can take place in the oceanic circulation and sea ice in the Bering Sea. Seasonal and interannual variations of sea ice in the Bering Sea have significant implications for shifts in marine species composition and ecosystem reorganization in the Bering Sea and adjacent oceanic regions [30,31]. Comprehensive investigations are required to gain a reliable knowledge of the aforementioned variations. In the present study, we employ the sea ice model coupled with a ROMS using a high resolution to investigate these multiple-scale variations in the sea ice and ocean circulation in the Bering Sea. This paper is organized as follows. Section 2 introduces data sources along with data inputs, Section 3 describes model configuration, simulation results are given in Section 4, and Section 5 presents the summary and discussion.

2. Data Sources

The following observational measurements are used in the present study: sea ice concentration, sea surface temperature (SST), eddy kinetic energy (EKE), sea surface wind, and air temperature. Sea ice concentration data were obtained from the 1990–2004 monthly SSM/I-measured (Special Sensor Microwave/Imager) microwave brightness temperature data, downloaded from https://nsidc.org/data/NSIDC-0051/versions/1, which is maintained by the National Snow and Ice Data Center (NSIDC) [32]. The data in the polar stereographic projection are provided with a grid cell size of 25 × 25 km. The grid size in the coordinates used is not exactly 25 km (see https://nsidc.org/data/polar-stereo/tools_geo_pixel.html) but slightly smaller than 25 km. The maximum difference is about 4% and the average difference is 2%. Therefore, on average, the sea-ice area and extent calculated from the SSM/I data is about 8% larger than the true sea-ice area and extent. The data were generated using the NASA Team algorithm, which uses a polarization ratio at 19 GHz and a gradient ratio of 19 V and 37 V. The average accuracy of NASA Team algorithm is within $\pm 5\%$ in winter in a compact ice pack, but it is hard to validate datasets at 0% and 100% because they are not designed to enable retrievals outside a sea ice concentration range of 0%–100% [33].

AVHRR-v2(AVHRR stands for “Advanced Very High Resolution Radiometer data”.) SST data is also used (https://www.ncei.noaa.gov/thredds/catalog/OisstBase/NetCDF/AVHRR/catalog.html). The data has a daily temporal resolution and a 0.25° grid resolution. The EKE was calculated from the AVISO(AVISO stands for “Archiving, Validation and Interpretation of Satellite Oceanographic data”.) altimetry-measured Sea Surface Height Anomaly (SSHA) (downloaded from http://www.aviso.oceanobs.com/duacs/), with a spatial resolution of $1/3°\times 1/3°$ and a 7 day temporal window over the period from January 1993 to December 2004.

Monthly sea surface wind and air temperatures from NCEP (National Center of Environmental Prediction) reanalysis dataset were used to analyze the physical mechanism. The 6 hourly (four times a day: 0000, 0600, 1200, and 1800 UTC) meteorological data from the same data source were used to force the numerical model. All the data were downloaded from https://rda.ucar.edu/datasets/ds093.0/ which has a spatial resolution of 0.3 degree.

3. Model

3.1. ROMS

ROMS solves rotating primitive equations with a realistic equation of state [34]. The model uses a generalized sigma-coordinate system in the vertical direction and a curvilinear grid in the horizontal plane. ROMS is a split–explicit, free-surface oceanic model, in which short time steps are used to advance the surface elevation and barotropic momentum equations, with a larger time step used for temperature, salinity, and baroclinic momentum. A third-order, upstream-biased advection operator allows the generation of steep gradients in the solution, enhancing the effective resolution of the solution for a given grid size when the explicit viscosity is small. The vertical mixing is parameterized using a K-profile parameterization (KPP) scheme [35].

3.2. Sea Ice Model and Processing

A sea ice model coupled with ROMS is implemented with sea ice dynamics a and thermodynamics module [29]. The sea ice dynamics are based on an elastic-viscous-plastic (EVP) rheology from Hunke and Dukowics [36] and Hunke [37]. The sea ice thermodynamics are based on Mellor and Kantha [38] and Häkkinen and Mellor [39], with two ice layers and a single snow layer used to solve the heat conduction equation. The atmospheric heat flux is applied to the ice surface, and the stresses caused by the wind and oceanic currents are incorporated into the sea ice model on the top and bottom of the sea ice, respectively. The salinity flux is imposed onto ROMS when the sea ice is formed or melted. Since the Boussinesq approximation is used in ROMS, which ignores density differences except terms multiplied by acceleration due to gravity, water mass is conserved. Therefore, the water mass change due to sea ice melting and freezing is ignored. The output variables include sea ice concentration (percentage per grid) and ice thickness (grid-cell mean). The monthly means of these variables are calculated based on 5 day averaged outputs by averaging 6 records of 5 day averages.

3.3. Model Configuration

The model domain, covering the whole Bering Sea by extending southward to cover the Aleutian Stream and northward to cover the southern portion of the Chukchi Sea, is plotted in Figure 1. The grids follow the zonal and meridional directions, with a maximum grid size of approximately 7 km on the southern end and 3 km on the northern end. West of Kamchatka Peninsula is considered land and is masked off since it has a neglectable influence on the Bering Sea. At the surface, the model is forced by the 6 hourly NCEP reanalysis product with a spatial resolution of 0.3 degrees [40]. The following meteorological parameters used for the forcing are sea surface wind at 10 meters, net heat flux (shortwave radiation and longwave radiation), cloud coverage (albedo), precipitation, and evaporation. The lateral fluxes along the open boundary are interpolated, as interpreted from the monthly SODA (Simple Oceanic Data Assimilation) product. The restoring data for the lateral open-boundary conditions are from the 1990–2004 monthly SODA global oceanic reanalysis product, with a horizontal resolution of $0.5\times 0.5$ degrees and 40 vertical levels [41,42], downloaded from http://iridl.ldeo.columbia.edu/SOURCES/.CARTON-GIESE/.SODA/.v2p0p2-4/. SODA data are applied to the model grid, including temperature, salinity, currents, and Sea-Surface Height (SSH). The solid boundary around islands and the mainland features no-normal and no-slip conditions implemented through a land-mask algorithm [43].

In addition to the meteorological flux at the sea surface, tidal forcing is applied to the simulation along the open boundaries. The tidal amplitudes and phases in both sea surface heights and the barotropic velocities of eight tidal constituents are obtained from a global inverse barotropic tidal model (TPX06) [44,45], which has a horizontal resolution of $0.25°$. The eight tidal constituents are $M_2$, $K_1$, $O_1$, $S_2$, $N_2$, $P_1$,$K_2$, $Q_1$, ordered by their amplitude in the region, and their definitions refer to the website of https://tidesandcurrents.noaa.gov/glossary. The barotropic transport from the TPX0.6 solution is adjusted based on the ROMS bathymetry because the bathymetry fields from the ROMS and TPX06 are different. To account for the 18.6 year cycle of astronomical tidal-generating potentials, nodal correction is applied to the sea surface height and the barotropic transport of the TPX06 [46,47]. These two procedures proved to be important to achieve an accurate barotropic tidal simulation [46].

The model starts from the state in December of 1990 from the SODA product produced by 1990 NCEP forcing and repeated for five years until the model reaches a quasi-steady state. Then, the model is run for 15 years (1990–2004). The five-day averaged model outputs are analyzed in the present study.

4. Results and Analysis

4.1. Seasonal Variation

As introduced in Section 1, the sea ice in the Bering Sea follows a seasonal variation. The Bering Sea is essentially ice free in July through September. Only a few remnants of sea ice are present in June, and sea ice starts to form in October. This can be clearly seen from the remote sensing data of the sea ice.

To calculate the monthly sea-ice area in the Bering Sea from 1990 to 2004, we multiply the sea-ice concentration with the grid cell area at each grid for each month, and then sum all the grid sea-ice areas. The time series of the 15 year averages of monthly sea-ice areas in the Bering Sea for each month are plotted in Figure 2a in a solid line. Dashed lines are the values plus/minus one standard deviation representing its interannual variation. The figure shows that the sea ice area reaches the maximum in March and the minimum in August and illustrates that the sea ice starts to form in late September and early October. The standard deviation (the difference between the dashed lines and solid lines in Figure 2) represents the interannual variation for each month. It can be clearly seen that the interannual variation varies from month to month. In winter and early spring (DJFMA), it is much larger than that in summer and early fall (JASO). ROMS accurately reproduces the seasonal variation in the sea ice coverage shown in the Figure 2b. However, ROMS produces less sea ice than the SSM/I observed during the winter peak because melting processes are delayed during early summer. The ROMS results exhibit a full month delay. Zero sea-ice areas are found in August in ROMS, while they appear as early as July in the observations. The largest differences seem to occur in May and June. Further, the freeze-up is delayed, as shown by the smaller sea-ice area in November and December in ROMS compared to the observations.

To check the spatial distribution of the sea ice coverage in the Bering Sea, the mean distributions of the sea ice concentration in March and August (averaged over the 1990–2004 SSM/I data) are plotted (Figure 3a,b). In March, the sea ice can be found covering most of the northern shelf of the Bering Sea. On the western Bering Sea, a narrow extension is presented along the eastern coast of Siberia and even northern part of the Kamchatka Peninsula. On the eastern Bering Sea, the whole western coast of the Alaska is covered by the sea ice but with a relatively lower concentration than in the western open ocean. In August, the sea ice in the Bering Sea completely disappears. It should be noted that, to maintain consistency, the land masks for both the model and remote sensing observational data (including sea ice, SST, and AVISO SSH data) are the same ones from the model.

The sea ice model coupled into ROMS can reproduce the above seasonal variation in the sea ice concentration in the Bering Sea. Figure 3c,d plots the model-simulated mean sea ice concentration distribution in March and August, also averaged over 1990–2004. The ice edge produced by the model is very close to that from the SSM/I observations in all directions. There are some differences in the concentration values between the model and observations. Figure 3e,f plots the differences between ROMS and SSMI obtained by ROMS minus SSMI. As can be seen, for example, there is a dipole in the difference in March, which tells us that the ROMS produces more sea ice along the western coast in the southern part and less sea ice along the Bering slope. The ROMS model also fails to reproduce high ice concentrations close to the coast towards the Kamchatka Peninsula. The reasons for these differences might be caused by the imperfect numerical models used. Numerical model systems, such ROMS and the sea ice module, are complicate systems, in which many factors could produce errors, such as the inaccuracy of the ROMS numerical schemes, which calculate the physical processes, such as turbulent mixing and sea-ice-ocean interactions. Some differences could be caused by resolution differences in the SSM/I data and the model configuration. To explore the reasons for such discrepancies requires a series of numerical sensitivity experiments, which is beyond the scope of the present paper.

SST data from satellite observations of the ice-free ocean surface are compared with the model solution. Figure 4a–d show the monthly anomalies from ROMS and satellite AVHRR SST data, respectively. The anomalies are calculated by using a 15 year monthly mean SST minus a 15 year yearly mean SST. Both the model and observations show opposite anomalies during spring and summer. In spring, along the Aleutian Islands and the Bering Slope, positive anomalies are seen for both the model and data, but negative anomalies are observed in summer. Along the Kamchatka coast, the anomalies are different from those along the Bering Slope. In spring, negative anomalies are present while positive anomalies occur in summer. Though the spatial distribution pattern of the observed SST is reproduced by the model, there is systematic discrepancy between the modelled results and satellite observational data. Such numerical errors could be caused by the uncertainties of numerical schemes implemented in the ROMS (e.g., air-sea interaction and upper layer mixing schemes), which need be improved in future.

The oceanic circulation in the Bering Sea also has a correspondingly large seasonal variation. Figure 5 displays the surface mean circulations in March and August averaged over 1990–2004. In March, the circulation is characterized by three jets in the Southern Bering Sea. A westward Aleutian stream enters the Bering Sea through gaps among the islands, a strong jet flows northwestward along the Bering Slope, and a current flows southward along the Kamchatka coast. Sea water spreads over the slope and onto the shelf at the north of the Bering slope, where the sea ice is presented during the wintertime. The poleward surface current is very weak across the Bering Strait during the winter. The circulation pattern shows a dramatic change, especially over the shelf north of the Bering Slope, where the sea ice is absent in August. It is also noted that poleward transport in the Bering Strait is observed much more strongly in August than in March.

The monthly poleward transport through the Bering Strait from ROMS is plotted in Figure 6, in which the solid line is the 15 year mean and dashed lines are values plus/minus one standard deviation value displaying the interannual variation. This demonstrates that the water exchange between the Pacific Ocean and the Arctic Ocean occurs mainly as the Pacific water flows into the Arctic Ocean, which is consistent with reports from the literature, such as the studies by Woodgate and Rebecca [48], Aksenov et al. [49], and Watanabe [50]. Figure 6 shows a strong seasonal variation of poleward transport. In spring and summer, the poleward transport peaks. In November through January, its magnitude is much smaller compared to summer. The poleward transport also carries strong interannual variation, as shown in the figure. The interannual variability (denoted by the standard deviation) is smallest for April through July and largest for November through February. It is 2–3 times more likely to have a poleward transport in February, which is as large as the one in May, than to have a poleward transport in May, which is as low as the result in February. From November through January, the poleward transport reversed to be southward from the standard deviations presented in Figure 6.

The seasonal variation in the sea ice coverage, oceanic circulation, and poleward transport can be explained by sea surface wind forcing and air temperature. Figure 7 shows the sea surface wind and air temperatures in March and August, averaged over 15 years from 1990 to 2004. In March, the strong northeasterly wind dominates, which carries cold air (a few degrees below 0 ℃) from the Arctic Ocean to the Bering Sea, causing sea ice formation in its northern part. The mean wind speed magnitude reaches about 10 m/s. This strong wind facilitates a southward spread of sea ice. In August, a mild southwesterly wind (weaker than 3 m/s) blows over the Bering Sea with air temperatures of about 5–10 ℃ [17]. This temperature helps melt the sea ice completely in the Bering Sea. It should be noted that a cyclonic center of the wind is clearly seen southeast of the Kamchatka Peninsula on the wind map in March, located near the Aleutian Low, but the cyclone disappears in August.

It can be understood that sea surface wind and air temperature changes are related to large- scale climate pattern changes (Aleutian Low and Siberian High) and drive the seasonal variation in sea ice and oceanic circulation [51,52].

4.2. Interannual Variation

The standard deviation in Figure 2 is calculated based on 15 year monthly means and represents the interannual variations in the sea-ice area. As can be seen from Figure 2, the interannual variation varies from month to month. To further demonstrate the interannual variation, Figure 8 shows the yearly means of the sea-ice area from both the model (ROMS) and the observations (SSM/I), for 15 years from 1990 to 2004. This figure clearly shows a decreasing trend in the sea-ice area during the 15 years from 1990 to 2004 for both the ROMS run and the satellite data, with about a 15–20 percentage decrease. The yearly time series also shows a 3–5 year oscillation. The decreasing trend in the sea-ice area within the studied 15 years could be related to a decadal oscillation in a Pacific domain scale [53]. It should be noted that there is relatively larger difference between the model and observation during the period from 1998 to 2002 than in other years. This difference could be caused by the uncertainties in the forcing for the model and the observational data.

As discussed in Section 4.1, seasonal variations in the sea ice and the surface oceanic circulation in the Bering Sea are partially driven by changes of the sea surface wind and air temperature. The latter is highly related to the Aleutian Low. What drives the interannual variation? To address this question, the correlations between several variables are determined (Figure 9, Figure 10 and Figure 11), including:

• The sea ice coverage anomaly and sea surface air temperature;
• The poleward transport anomaly across the Bering Strait and the wind intensity anomaly;
• The sea ice coverage anomaly and Siberia-Aleutian Index.

The time series of the domain-averaged air temperature from the NCEP and the sea ice coverage with the seasonal cycle removed are plotted in Figure 9, and the correlation between air temperature and sea ice coverage is −0.66 at the significance level of 95%. This shows the sea ice variation is significantly affected and correlated with more (less) sea ice, which associates with lower (higher) air temperature.

The domain-averaged poleward transport is defined as follows:

$$T_y\left(y,t\right)={{\displaystyle \int }}_0^{Lx}{{\displaystyle \int }}_0^{H}v\left(x,y,z,t\right)dxdz$$

$$T_{total}\left(t\right)=\frac{1}{Ly}{{\displaystyle \int }}_0^{Ly}{T}_y\left(y,t\right)dy$$

where $Lx$ and $Ly$ are the zonal and meridional sizes of the numerical model domain, respectively. $\mathrm{H}$ is the water depth.

Figure 10 displays the time series of the domain-averaged poleward transport variation and poleward wind speed variation. The correlation reaches 0.71 at a significance level of 95%. This high correlation demonstrates that wind forcing is the main driving factor, which is consistent with Yang and Dai [54]. They applied a coupled model to investigate the effect of the ocean surface winds on meridional transport and found that the absence of surface wind causes a remarkable decrease in poleward oceanic transport. The interannual variation in Figure 8 can also be observed in the poleward mass transport across the Bering Strait, shown in Figure 10.

The Siberia-Aleutian Index represents the difference between the mean winter (DJFM) normalized 700 hPa anomalies in two regions, Siberia (55$°$N–70$°$N, 90$°$E–150$°$E) and Alaska/Yukon (60$°$N–70$°$N, 130$°$W–160$°$W). The base period for the index normalization is 1961–2000. The 700 hPa height data were obtained from the NCEP Reanalysis Project (https://rda.ucar.edu/datasets/ds093.0/). The index represents the local wind intensity in the area and is used to describe the variation of the local climate system. Figure 11 plots the time series of the Siberia-Aleutian index and the sea ice coverage variation. Their correlation is 0.68 at a significance level of 95%. The Siberia-Aleutian index represents a large-scale wind pattern that drives sea surface current and ocean mixing. The oceanic current and the upper ocean mixing could cause the SST variation, which directly affects the sea ice variation [55].

4.3. Intraseasonal (Eddy) Variation

It has been reported in the literature that there are strong eddy activities along the Bering Slope due to the instability of the strong slope current [1,2,3,4]. Five-day averaged model data are used to study the eddy or the intraseasonal variation.

A snapshot of surface vorticity is plotted in Figure 12, showing a well-developed eddy street, which is a repeating pattern of swirling vortices induced by the unsteady separation of the flow along the Bering Slope. The eddies are advected westward by the Bering Slope Current along the slope. Figure 13 plots the EKE calculated from both the ROMS and the AVISO altimetry-measured SSHA data, both of which show an EKE band along the Aleutian Islands, the Bering Slope, the Kamchatka coast, and the Bering Strait. However, the magnitude of the EKE differs between the ROMS and AVISO measurements. This discrepancy is caused by many factors, including the inaccuracy of the model’s numerical schemes, simulating physical processes in the ocean, such as uncertainties in the mixing parameterization, errors in the satellite measurements, the resolution between the ROMS run, and the satellite observations.

Eddy momentum transport across the Bering Slope is plotted in Figure 14b. The seasonal signals larger than 90 days have been high pass filtered. To calculate the eddy transport across the Bering Slope, the coordinates were rotated $\left(90-\mathsf{\theta }\right)$ anticlockwise, where $\mathsf{\theta }$ is the angle between the Bering Slope and the negative direction of the X axis. The velocity components $\mathrm{u}$ and $\mathrm{v}$ are transformed to ${\mathrm{u}}_{\mathrm{r}}$ and ${\mathrm{v}}_{\mathrm{r}}$, where ${\mathrm{u}}_{\mathrm{r}}$ is perpendicular to, and ${\mathrm{v}}_{\mathrm{r}}$ is along, the Bering slope. The magenta straight line approximately follows the 1500 m isobath along the Bering Slope (as shown in Figure 14a). It can be seen that eddy momentum transport is significantly affected by the Bering Slope and mainly located over the Bering Slope. It is also noted that the eddy momentum transports are also large along the Aleutian Islands, the Bering Slope, the Kamchatka coast, and the Bering Strait, which is consistent with the EKE distribution.

4.4. Polynya Evolution

Polynyas are either oceanic areas of open water surrounded by sea ice or areas covered with new or young ice [56]. The St. Lawrence Island polynya is a common phenomenon in the Bering Sea and occurs at least once annually during the winter season [57,58].

As an example, Figure 15 presents the time series of the 5 day averaged sea ice concentration variation around the St. Lawrence Island in April and May, 1996. Polynya evolution can be clearly seen in Figure 15. On 5 April, the polynya begins to form along the southern coastline of the St. Lawrence Island, and increases southwestward when it reaches the maximum extent on 25 April. In May, the polynya forms in the northern coastline of the St. Lawrence Island.

The St. Lawrence Island polynya is formed by the sustained northerly winds in the Bering Sea [57,59], and the St. Lawrence Island polynya is classified as a latent heat polynya with dense water formation south of St. Lawrence Island. Dense water formed during the polynya ice production is the key factor. Therefore, polynyas generate and maintain open-ocean deep convection, which influences thermohaline circulation [57,60,61,62].

5. Summary and Discussion

In this study, the multiple-scale variations of the sea ice area and oceanic circulation in the Bering Sea (intra-seasonal (eddy), seasonal, and interannual) are investigated. Both 15 year (from 1990–2004) satellite remote sensing data and numerical modeling results are employed in the study. This model can reproduce the seasonal and interannual variation in the sea ice area compared with the satellite remote sensing data, including the SST, altimetry-measured SSH, and sea ice. The spatial distribution of the SST and its seasonal variation is exposed by remote sensing data, which is reproduced by a ROMS simulation. The seasonal variation in the sea ice coverage composes the primary variation in the Bering Sea, while the interannual variation is also present, which is revealed by both the model and satellite remote sensing data. Strong intraseasonal variation or events (polynya) can be found around islands via the high-resolution numerical modeling results. The oceanic circulation also has seasonality. A strong vortex street, which separates the shallow part from the deep part in the Bering Sea, is present along the Bering Slope. This sheet can be seen from both the AVISO data and the numerical modeling results.

To explore the mechanisms driving the multiple scale variations in the Bering Sea, the correlations between sea ice coverage and different dynamic parameters, such as the air temperature, SST, sea surface wind, and basin-scale Siberia-Aleutian index are calculated based on both numerical modelling results and observational data. It is revealed that these above variables are important in controlling sea-ice variations in the Bering Sea.

However, there are deviations between the numerical model results from the observational data. These differences could be caused by either uncertainty in the numerical schemes used in the numerical model or observational errors. For examples, there are many empirical parameterizations that need to be improved in both the oceanic numerical model and the sea ice model. The technique used for the remote satellite sensing of sea ice coverage and the algorithm in the calculation of the sea ice need to be improved furthermore. The differences in the model–observation comparison could result from the different resolutions of the two data. Therefore, more effort is needed to improve both the sea-ice ocean coupled model and the observational technique.