Characteristics of Atmospheric Circulation Associated with Variability of Sea Ice in the Arctic

Abstract

The paper investigates the role of atmospheric circulation in the surface layer in forming the Arctic ice structure. For the analysis, the empirical orthogonal function (EOF) method of decomposition of the surface wind field is used, and the reaction of ice to changes in the principal components of leading EOF modes is investigated using statistical methods. Analyzing the rate of ice change in the Arctic associated with the Arctic ocean oscillation mode, we concluded that this mode’s variability leads to the formation of a seesaw in the ice field between two regions. From the one side, it is the region of the central deep-water part of the Arctic, including the East Siberian Sea, and from the other side, it is all other marginal seas. The second (“dipole”) mode is most associated with an increase/decrease in the ice thickness at the Arctic exit through the Fram Strait, as well as the formation of the so-called “ice factory” in the coastal region of the Beaufort Sea in the positive phase of this mode. There is also a significant relationship between the variability of third mode and the arrival of Atlantic waters with a high heat content into the Arctic through the Barents opening, which creates preconditions for ice formation in this region.

1. Introduction

The variability of the Arctic sea ice distribution is mostly determined by atmospheric circulation since it drives the sea ice advection and forms inhomogeneities in the ice movement (dynamic factor). It is also responsible for transporting heat, cloudiness, and humidity by air currents, determining the balance of heat flows at the ice surface (thermodynamic factor) [1,2]. Feedback is determined by the presence of ice-free areas, the change in the albedo of the underlying surface, and the change in the surface air temperature and humidity. However, we should emphasize that the atmospheric factors mainly work in addition to the ocean factors such as ocean thermal anomalies, surface currents, advective heat transport, and the structure of the upper thermohaline.

The dominant mode of atmospheric circulation in the Arctic and adjacent mid-latitudes is the Arctic oscillation (AO, see the list of abbreviations at the end of the article) [3]. Since the 90s, this has become even more evident. It was shown that the interannual variability of the northern winter stratospheric flow in 1964–1993 was closely linked to large-scale circulation anomalies in the middle troposphere [4]. Using the EOF method of decomposition (Empirical Orthogonal Functions) between the 500 and 50-hPa geopotential heights, a 500-hPa AO structure was produced, including anomalies in eastern Siberia. Based on numerical modeling in [5], the relationship between tropospheric and stratospheric circulation and ocean surface temperature was established.

A noticeable decrease in the volume and area of sea ice [6] was associated, on the one hand, with the dynamic features of the atmospheric circulation [7], on the other, with the influence of warm air masses [8]. Thanks to [9,10], the connection with the general structure of the atmospheric circulation, represented to a large extent by AO, became more apparent.

Even though the movement of ice and its piling (pressure ridges), caused by the wind’s action, do not directly change its volume, nevertheless, the consequences of such movement can be significant. For example, an increase in ice transport along with the Transpolar drift leads to its more rapid melting at lower latitudes, and a weakening of this drift, on the contrary, leads to the accumulation of multi-year ice in the Arctic.

Another example is the lack of imported ice in the Barents Sea in 2010–2016. Ice deficit leads to a shortage of surface freshwater formed during thawing, weakening of stratification, and, as a result, to an increase in vertical mixing, rise of warm and saline water to the surface, and weakening of the ice formation process in the region [11]. Thus, the vertical structure of the Barents Sea waters is transformed from a cold stratified structure of the Arctic type into a warm well-mixed structure characteristic of the North Atlantic.

The action of the thermodynamic factor directly affects the volume of ice. A positive or negative balance of heat fluxes on the upper, lower, or lateral ice surfaces leads to ice melting or seawater freezing. Snowfall increases the thickness of the snow cover, and later, as a result of deformation and fusion of snowflakes, the snow transforms into ice mass.

The thermodynamic factor also becomes important in the final phase of dynamic ice movements. If the thickness of the ice changes due to advection, then, as a rule, this leads to heat fluxes imbalance. Ice that has become too thick due to advection begins to melt, and ice that has become too thin begins to freeze. As a result, we can expect from a long-term perspective that the dynamic and thermodynamic factors will act in contrary directions when the balance is reached.

Using statistical methods and the EOF decomposition method in several previous works (for example, Ref. [12,13,14,15]), they showed a direct relationship between atmospheric circulation modes and the nature of ice distribution and its dynamics in the Arctic. Based on observations from 1953–1992, Ref. [14] proposed a mechanism for the relationship between the anomalies in ice concentration and sea level pressure, ensuring cyclic repetition every ten years. The increasing trend of the North Atlantic Oscillation (NAO), according to [10], caused a significant reduction in ice in the Arctic Ocean.

The existence of negative feedback for the time scales of several weeks was shown in [16] using numerical simulation. The initial reaction to the anomaly in ice concentration and surface temperature associated with the NAO’s positive polarity consists of an opposite development of the geopotential anomalies in the lower and upper troposphere, which form near the disturbance area. Thus, the negative feedback process begins with a localized baroclinic response that reaches peak intensity in 5–10 days and persists for 2–3 weeks, gradually growing into barotropic and increasing in time and space. The equilibrium state is reached after 2–2.5 months. The atmosphere develops a larger-scale, equivalent barotropic response that resembles the NAO’s negative polarity, and this pattern is maintained primarily by nonlinear transient eddy fluxes of vorticity, related in part to changes in tropospheric Rossby wave breaking [17]. Statistically, such a development of the atmospheric response was also confirmed by [18]. Under global warming conditions, the noted development is likely to be similar but will be modified due to the reduction of the ice field [19].

In this article, we present some results of an analysis of the sea ice thickness dynamics, the rate of its change due to ice movement and thermodynamics, and how this characteristic develops depending on the values of the main atmospheric circulation indexes, the most important of which is the Arctic Oscillation (AO) Index.

2. Model, Data and Numerical Experiment

To estimate the characteristics of Arctic sea ice dynamics and thermodynamics, we used the simulation results obtained using a coupled ocean-ice model of the Arctic and the North Atlantic–SibCIOM, Siberian Coupled Ice-Ocean Model. Previously, it was referred to as the ICMMG model (Institute of Computational Mathematics and Mathematical Geophysics, Novosibirsk, Russia) by the institute’s name, which developed it. The model coordinate grid was constructed using a three-polar system and had 1/2${}^{\circ }$ resolution in the North Atlantic and 10–26 km resolution in the polar regions. The model has 38 vertical levels. The upper 400 m layer contains 18 levels with a resolution of 5 m at the surface to 50 m at 400 m. The upper 1000 m layer includes 24 levels.

The model used PHC data [20] for initialization. Sea ice is described using the CICE-3 model and includes five ice categories and one category of snow. The initial ice was 2 m thick and had a 100% concentration where the starting surface temperature of the ocean was below zero. After that, the coupled model run for the period from January 1948 to December 1956. All subsequent experiments used the resulting ice field to set the initial distribution. The IBCAO [21] bathymetry was a basis for constructing the model bathymetry. The river discharge of the region’s main rivers was taken from the data of RivDis-1.0 [22] and annually repeated the averaged seasonal variation. Besides, we were applying an increase or decrease in river flow to offset the imbalance between evaporation and precipitation. The Bering Strait discharge was set equal to 0.8 Sv and from 1998 to 2008, following [23]. At the southern boundary (20S), a uniform discharge was set, compensating for all external tributaries.

Most of the results used in our analysis we obtained during a numerical experiment of restoring the ocean and sea ice dynamics from January 1948 to December 2019. Atmospheric forcing was obtained during the reanalysis of NCEP/NCAR [24] and includes the temperature of surface water (at $\sigma =$ 995 db), humidity, sea level pressure, precipitation rate, descending long-wave and short-wave radiation, as well as wind speed in the surface layer.

The previous studies proved the possibility of the use of this model. Based on numerical experiments conducted using atmospheric forcing of the modern reanalysis of NCEP/NCAR, CORE-II [25], we performed an analysis of recent climate changes in the Arctic Ocean [26,27,28]. We also considered the issues of: (a) model reproduction of observed ice thickness and compactness variations, (b) the influence of atmospheric circulation variations on ice drift, and the circulation of the upper ocean layer, (c) changes in the thermohaline structure of the Arctic Ocean waters, and the thermodynamics of the ice cover [29,30] also showed the influence of Pacific and Atlantic waters on the distribution and Arctic ice thickness.

3. The Main Modes of EOF Decomposition and Their Characteristic Features

3.1. EOF Decomposition Method

To study climatic fields, we will use an empirical orthogonal function (EOF) analysis. Specifically, each eigenvector of the spatial covariance matrix of the data is the corresponding EOF (or EOF mode), and the time series obtained by projection of the data onto that EOF is its principal component. The method allows us to reduce the dimension of the analyzed system and find a relatively small number of independent variables that represent the essential information about the variability of the source data. The method also allows one to evaluate the relative importance of each structure.

First, the available data must be organized in the form of a matrix $Z_{ij}$. The columns of this matrix are a “snapshots” of the data for all stations $i=1,\dots ,N$ (space grid points) at the corresponding time moment $t_j$, and the rows are the time series of observations $j=1,\dots ,M$ for some station $(x_i,y_i)$ (grid point). EOF analysis uses a set of orthogonal functions to represent data as follows:

$$Z_{ij}=Z(x_i,y_i,t_j)\approx \sum _{k=1}^K{P}_k\left(t_j\right)E_k(x_i,y_i),$$

where $E_k$—empirical orthogonal function, $P_k$—its corresponding principal component [31]. Functions are calculated sequentially and each k-th of them is selected orthogonal to the previous $(k-1)$ functions so as to minimize the residual difference when expanding in $(k-1)$ functions. Since we analyze gridded data, each of $E_k$ represents a spatially distributed function—a k-th mode.

Using monthly averaged of NCEP/NCAR reanalyses data for each January, February, and March from 1948 to 2019, we define a new index as the leading principal component $P_1\left(t\right)$ of sea level pressure over the northern hemisphere starting from 20 N:

$$Z_{ij}=\left(p_i\left(t_j\right)\right),i=1,\dots ,N,\mathrm{such} \mathrm{that} y_i\ge 20\mathrm{N}.$$

Figure 1a shows an associated EOF and features the well-known Arctic Oscillation structure.

A characteristic feature of the distinguished leading mode is the anti-phase behavior of polar latitudes (>60 N) and moderate (>20 N and <60 N), leading to an increase or decrease in the polar vortex. The structure of this mode in the North Atlantic region is such that the pressure difference between the Icelandic minimum and the Azores maximum also largely depends on this mode, which ensures the relationship between AO and the North Atlantic Oscillation (NAO) [3,32].

Figure 1b presents time changes in the corresponding principal component in comparison with the AO index, according to the data presented on the website http://www.cpc.ncep.gov (12 April 2020). Both time-series were normalized so that their average and r.m.s. values are equal accordingly to 0 and 1. As a result, it turned out that our index correlates with the AO index with a coefficient of 0.9852. To make sure that such a high value is not a consequence of the general trend, we subtracted the quadratic trends of both data sets. The correlation coefficient for detrended data slightly decreased to 0.9624 but is still high. This high value allows us to assert that the method we developed for determining the oscillation index corresponds to the generally used AO definition, and the visual similarity of the eigenfunction with the Arctic Oscillation structure turned out to be quite natural.

Speaking about the correlation coefficient r, one must keep in mind the significance of this or that value. If number t obtained as a result of the conversion

$$t=\frac{r}{\sqrt{1-r^2}}\sqrt{M-2},$$

where r is the correlation coefficient, M is the sample length, then the sample distribution of t is a Student’s t-distribution with $(M-2)$ degrees of freedom. For the analyzed time-series, the sample length is $M=$ 864 (72 years by 12 months). If we take the absence of correlation as the null hypothesis, i.e., $r=0$, then with a confidence of 95%, it can be rejected if $\left|t\right|>t_{cr}(M-2)$. For such a large M, the value of $t_{cr}(M-2)$ is 1.96. It means that the correlation coefficient is significant if $\left|r\right|>$0.07, which is a rather soft requirement. In our analysis, we will use a sliding average with a 5-year averaging period. Assuming that we thereby reduce the sample length by 60 times (5 years by 12 months), we get the new sample length $M^{*}=$ 14, with $t_{cr}(M^{*}-2)=$ 2.1788 and a significance threshold for correlation coefficient rises to the level of $\left|r\right|>$0.53. However, the use of a moving average does not mean an unambiguous decrease in the degrees of freedom corresponding to the averaging period, since we still take into account all 864 values of the series and do not reduce their number down to 14 values. The true value of the criterion should be somewhere between 0.07 and 0.53. We assume that the critical value must be about 0.5, which is close enough to the criterion for a sample of 14 values. Accepting this, we will further say that the value of the correlation coefficient is significant if the condition $\left|r\right|>$0.5 is correct, which we assume to mean that with a confidence of 95% the hypothesis that there is no linear relationship can be rejected.

3.2. EOF Decomposition of the Surface Wind

When analyzing the variability of the ice field, the sea level pressure does not represent that part of atmospheric forcing that directly affects this variability. The most of sea ice changes are the results of thermodynamic processes (water freezing, ice melting), ice dynamics (convergence, divergence, ice ridging), precipitation (mainly in the form of snowing over the ice surface). In the previous work [28], we examined thermodynamic processes in detail and concluded that the most significant changes due to climatic trends in the ice field occur due to changes caused by long-wave radiation fluxes. We omit the details which can be found in [28], proving that with an increase in surface temperature, air humidity near the ocean surface increases, and as a result, a larger part of the long-wave radiation emitted from the surface is reflected by moist air, creating a semblance of the greenhouse effect.

To study the dynamic component of the ice field variability, we will consider the wind velocity vector in the atmospheric surface layer from the NCEP/NCAR reanalysis data, more precisely its two horizontal components: zonal u and meridional v at the level of 995 hP of air pressure (from now on we will refer to it as to a surface wind) for the period from 1948 to 2019. Besides, we will consider only the points of the Arctic Ocean and the adjoining seas of the Atlantic (Greenland and Norwegian), i.e., the part of the Atlantic, bounded from the south by the Danish Strait and by the line connecting Iceland, Farrero-Shetland Islands, and the Scandinavian Peninsula. Thus, the state vector in this case is

$$Z_{ij}=\left(u_i\left(t_j\right),v_i\left(t_j\right)\right),i=1,\dots ,N^{*},\mathrm{such} \mathrm{that} (x_i,y_i)\in \mathrm{\Omega },$$

where $u_i\left(t_j\right)$ and $v_i\left(t_j\right)$ are the surface wind velocity components at i-th point with coordinates $(x_i,y_i)\in \mathrm{\Omega }$ at time moment $t_j$, subscript $i=1,\dots ,N^{*}$ enumerates all NCEP/NCAR reanalysis data nodes inside area $\mathrm{\Omega }$, described above. In order to reduce the size of the matrix of states $Z_{ij}$, we use a monthly averaging over time, which, unfortunately, leads to filtering out daily and synoptic variability.

Figure 2a shows the surface wind field averaged over the indicated period and its vorticity. Basically, it demonstrates two cells of negative and positive values separated along the line of the Lomonosov Ridge. Negative (anticyclonic) rotation centers in the Canadian basin, and its extremum falls on the Beaufort Sea region. The area of positive (cyclonic) vorticity is in the region of the Amundsen Basin, the Barents, Kara and Laptev seas, and its center of the action, apparently, is located in the northern Atlantic seas, and is possibly associated with the entire North Atlantic [33], which is outside the scope of our analysis. Characterizing this field, we can note that the prevailing wind promotes the movement of the Transpolar Drift, the removal of ice from the marginal seas of the Eurasian shelf, as well as the inflow of ice into the Barents Sea through its northern border. The latter contributes to the replenishment of surface freshwaters in this region and enhances the stability of stratification.

3.2.1. First Mode

The first eigenfunction (or the first EOF mode) is presented in Figure 2b, and its principal component, smoothed using 5-year sliding averaging and normalized after smoothing (that is, reduced to the form when the mean is zero and the standard deviation is unity), is presented in Figure 3a. The changes associated with this mode explain about 26% of the surface wind variability, and taking into account the sampling error, this mode explains, according to [34], 23–28.2%. The variability of this mode by approximately 21% is seasonal (i.e., the ratio of the standard deviation of its averaged seasonal variability is approximately 21% of the total standard deviation of this value). This percentage does not demonstrate clear seasonality. However, it follows from the analysis that this mode is more often in the positive phase in November–April, when the amount of ice is maximum, and in the negative phase in June–September, when ice is minimal.

In its structure, the eigenfunction of the first EOF mode resembles the distribution characteristic of the Arctic Ocean oscillation (AOO) [35], which characterizes two circulating regimes in the ocean—cyclonic and anticyclonic. From Figure 2b, one can see that the positive phase of this eigenfunction is associated with an increase in the anticyclonic regime in the Canadian basin, and in the negative phase, with its weakening and concentration in the Beaufort Sea region. This fact is in qualitative agreement with the results of [35]. The correlation coefficient of the principal component of the mode with the AOO index [36] is 66% (see

Table 1. Correlation (%) of the principal components of the surface wind EOF decomposition with the well-known circulation indexes: AO, NAO, AMO, AOO, PNA, PDO, taken from the site https://www.whoi.edu/page.do?pid=66578 (*—a fraction of the wind field variability explained by this mode, **—a fraction of the seasonal variability of the mode in its total variability).

Mode	% *	% **	AO	NAO	AMO	PNA	PDO	AOO
1	25.65	21.32	−41	−33	6	46	4	66
			(−73)	(−50)	(9)	(51)	(−4)	(74)
2	20.53	7.69	6	−31	75	46	−20	46
			(−29)	(−58)	(26)	(30)	(−40)	(20)
3	11.94	17.25	−77	−64	32	18	−32	34
			(−80)	(−57)	(30)	(54)	(−4)	(44)
4	8.29	10.19	−38	−55	42	−19	−53	−9
			(−22)	(−42)	(18)	(−6)	(−23)	(−25)
5	4.45	5.20	8	−11	56	60	12	48
			(−32)	(−32)	(20)	(46)	(−7)	(30)

), and after elimination of quadratic trends, it is even higher—74%. Therefore, hereafter, this mode will be called “oceanic”.

The relationship of this mode with the AO index at first glance seems insignificant (Figure 3a), and the correlation coefficient is only −41%. However, this is largely due to the presence of a temporary trend in the AO index. It is almost absent in the principal component. After eliminating the quadratic trend, the correlation coefficient for the residuals turned out to be −73%, which indicates the existence of a certain linear relationship between the two time-series. Among the other indexes presented in Table 1, the most consistent with this principal component is the Pacific North American Oscillation (PNA) index, which correlates with 46%. Let us assume that the AO and PNA indexes are independent (for example, Ref. [37] shows that the correlation between the Pacific and Atlantic centers is rather weak). Next, we make a linear combination in the form

$$\widehat{P_1}=\alpha P_A+\beta P_P,$$

where $P_A$ and $P_P$ are the normalized time-series of the AO and PNA indexes. Now we can pose the problem of finding such values of the coefficients $\alpha$ and $\beta$ that provide the minimum of the functional

$$J(\alpha ,\beta )=\sum _{j=1}^M{\left(P_1-\widehat{P_1}\right)}^2.$$

The solution to this problem gives the dependence of the coefficients $\alpha$ and $\beta$ on the following average values of time-series

$$\alpha =\frac{\overline{P_1{P}_A}·\overline{P_P^2}-\overline{P_1{P}_P}·\overline{P_A{P}_P}}{\overline{P_A^2}·\overline{P_P^2}-{\overline{P_A{P}_P}}^2},$$

$$\beta =\frac{\overline{P_1{P}_P}·\overline{P_A^2}-\overline{P_1{P}_A}·\overline{P_A{P}_P}}{\overline{P_A^2}·\overline{P_P^2}-{\overline{P_A{P}_P}}^2},$$

where the overline sign means averaging over a time-series. If we take 5-year sliding averages of the original series as $P_1$, $P_A$ and $P_P$, then the ratio between the optimal values of $\alpha$ and $\beta$ will be (−0.61):1, while the correlation coefficient between series $P_1$ and $\widehat{P_1}$ will be equal to 69% (compare with −41% and 46% for $P_A$ and $P_P$). Such a high correlation refers only to interannual variability and does not mean that these series behave synchronously on smaller time scales since we applied a sliding 5-year averaging. Therefore, for example, on the scale of seasonal changes, these series may not correlate at all. If, in addition to averaging, we apply the subtraction of quadratic trends, then the ratio of the optimal values of the coefficients will become (−1.9):1 (that is, the contribution of AO will increase by more than three times), and the correlation will increase slightly to the level of 76%, although for AO it was already quite high −73%. From this, we can conclude that the contribution of the AO in the first principal component of our decomposition is approximately two times greater than the contribution of the Pacific center of the action. However, when we consider the long-term trend (up to 70 years), the situation is the opposite—the role of the Pacific Ocean is more noticeable than the role of the AO. Since we stated before that the first principal component of our expansion coincides (up to 74% correlation) with AOO, then, speaking of the first principal component, we can also bear in mind AOO. However, a direct comparison of the resulting series $\widehat{P_1}$ with the AOO index gives only a 40% correlation.

AO is closely related to the North Atlantic Oscillation (NAO), for example, Ref. [33] believes that the AO is a global extension of the NAO. The correlation of the NAO index with the first principal component was lower than with AO, at the level of −33%. However, if we oppose to the Pacific center of the action, not the Arctic (AO), but the North Atlantic (NAO) center, and in the previous argument instead of the $P_A$ series we take $P_N$, which represents the normalized time-series of the NAO index, then the optimal ratio of coefficients for 5-year averaging after exclusion of a quadratic trend will be (−0.63):1. That is, the role of the Pacific Ocean is about one and a half times more significant than the role of the North Atlantic. Moreover, in the presence of trends, the correlation with $P_1$ is 61%, and without these trends, it is 65%. The correlation with the AOO time series is 45%, which is slightly higher than in the previous case.

A possible role of this mode in the formation of ice is to enhance or weaken the anticyclonic circulation in the Beaufort Sea region, i.e., to enhance or weaken the convergence of the Ekman transport (pumping) in this basin, which contributes to or counteracts ice consolidation.

3.2.2. Second Mode

The second EOF decomposition mode is shown in Figure 2c and Figure 3b. Its structure represents the distribution characteristic of the Arctic dipole [38,39], which is characterized by the pressure seesaw between North America (Greenland) and Eurasia (Taimyr). Thereafter, we will call this mode “dipole”. Figure 2c shows that the eigenfunction in its positive phase is associated with an increase in the Transpolar Drift. It contributes to the export of ice from the Arctic through the Fram Strait to the North Atlantic [38]. We can expect that this factor will have a significant impact on the dynamic component of changes in the ice field in the Arctic. The development of such a circulation structure is associated with climate change in the 2000s and represents northeastward shifts of the AO/NAO centers of action [40]. However, the correlation coefficient of the principal component of this mode with the AO index is only 6% (see Table 1), and after eliminating quadratic trends, it is slightly higher—−29% (anticorrelation). The time variation of the principal component does indicate a long period of the positive phase of the mode observed in the 2000s.

Among the indexes presented in Table 1, the most consistent with this principal component is the Atlantic Multi-decadal Oscillation (AMO) index, with a correlation of 75%. However, such a high indicator is a consequence of general trends, which can be either related to each other or merely coinciding. In favor of the former, for example, an analysis of observations and simulation data [41], confirming the existence of a connection between the trend of the AMO index and the Arctic Sea Ice Decline, proves. In the absence of trends, the correlation becomes as weak as 26% (Figure 3b).

In the absence of a long-term trend, the indexes closest to the second principal component are NAO and PDO (Pacific Decadal Oscillation). If we assume that these indexes are independent and, as in the previous consideration, make a linear combination of their normalizations with the coefficients $\alpha$ and $\beta$, then the ratio between the optimal values of $\alpha$ and $\beta$ will be (−4.8):1, with the correlation coefficient between the series $P_2$ and $\widehat{P_2}$ equal to 56%, which is slightly worse than if we take only the NAO index—−58%. This trend towards an increase in the values of the principal component coincides with the AMO index trend in the 2000s. In the absence of the trend, the correlation with the NAO index becomes more significant. Although small, it is significant.

3.2.3. Third Mode

Figure 2d and Figure 3c represent the third eigenfunction and its principal component smoothed using 5-year sliding averaging. The changes associated with this mode explain about 12% of the surface wind variability, of which approximately 17% are seasonal changes. From the analysis of its seasonality, it follows that this mode is more often in the positive phase in May–August when the ice melts intensively, and in the negative phase in December–March, when it freezes.

In its structure, the eigenfunction of the third EOF mode resembles the distribution characteristic of the NAO [33]. Unlike the AO structure, its action is mainly localized in the North Atlantic seas and further in the Arctic along the directions of Atlantic water flow through the Fram Strait and in the Barents Sea. In its positive phase, the anticyclonic vorticity of the wind increases in this region, creating a convergence zone here in the Ekman layer and counteracting the influx of Atlantic waters into the Arctic on the eastern flank and the outflow of Arctic waters on the western flank. Therefore, this mode will be called the “Atlantic” in the future. The maximum values of the correlation coefficient of this mode principal component are achieved with the index AO—−77%, and NAO—−64% (see Table 1), and after elimination of quadratic trends, respectively, equal to −80% and −57%. Among the Pacific indexes, after the elimination of trends, the PNA index stands out, with a correlation of 54%.

If we make a linear combination of normalized indexes AO and PNA

$$\widehat{P_3}=\alpha P_A+\beta P_P,$$

then the ratio between the optimal values of $\alpha$ and $\beta$ will be (−2.2):1, while the correlation coefficient between the series $P_3$ and $\widehat{P_3}$ will be 83% (in comparison with the values of −77% and 18% for $P_A$ and $P_P$). That is, the role of the North Atlantic in this mode is more than twice as significant as the role of the Pacific Ocean.

The manifestation of this surface wind mode is possible in enhancing the heat fluxes that additionally come along with the Atlantic waters in the negative phase while enhancing the removal of Arctic surface waters. However, the response of the ice field to these flows can be delayed in time, since for the Atlantic layer to significantly warm, the source must act for a long time, and besides, the thermal signal itself does not propagate instantly.

4. Results of Numerical Simulation

Analyzing the results of a numerical experiment, we will focus on such a characteristic as the growth/decrease rate of ice while considering the dynamic and thermodynamic growth both separately and in total. Since the type of atmospheric circulation determines the structure of the averaged surface wind, the dynamic component of ice growth is mostly determined by its mechanical transport under the influence of wind friction and ocean currents, as well as the divergence/convergence of the resulting ice velocity (Figure 4a). As a result of the emerging balance, thermodynamic growth, or melting (Figure 4b) should eventually compensate for the dynamical transport, so it is essential to consider this component as well. In general, Figure 4, obtained by numerical simulation, demonstrates that ice is formed in the central deep-sea regions of the Arctic and off the coast of the Laptev and Kara seas. Under the influence of a dynamic factor, it spreads from these regions toward the open Arctic and further towards the Greenland and Barents seas and towards the coastal part of Alaska and Chukotka and the Bering Strait, where its predominant melting takes place. This result is in good agreement with the available estimates [42,43].

We assume that the structure of atmospheric circulation is closely related to the leading EOF modes of surface wind decomposition considered earlier. In order to identify the characteristic features of the ice dynamics for a particular type of atmospheric circulation, we consider the averaged values of the ice growth rate for positive values of each of the series of principal components and separately for negative ones. Averaging will be performed, taking into account the values of the modulus of the principal components as a weight coefficient. Let $\widetilde{P_k}$ be the normalization of the principal component $P_k$ of the k-th mode of the EOF expansion after applying the sliding average, then we set the average ice growth rate in the positive phase $\overline{R_k^{+}}$ and in the negative phase $\overline{R_k^{-}}$ at each of the points $(x,y)$ equal to

$$\overline{R_k^{+}}(x,y)=\frac{{\sum }_{j=1}^M{W}_k^{+}\left(t_j\right)R(x,y,t_j)}{{\sum }_{j=1}^M{W}_k^{+}\left(t_j\right)},$$

$$\overline{R_k^{-}}(x,y)=\frac{{\sum }_{j=1}^M{W}_k^{-}\left(t_j\right)R(x,y,t_j)}{{\sum }_{j=1}^M{W}_k^{-}\left(t_j\right)},$$

where

$$W_k^{+}=\left\{\begin{array}{cc}\widetilde{P_k},\hfill & \widetilde{P_k}>0\hfill \\ 0,\hfill & \widetilde{P_k}\le 0\hfill \end{array}\right.$$

$$W_k^{-}=\left\{\begin{array}{cc}0,\hfill & \widetilde{P_k}\ge 0\hfill \\ \left|\widetilde{P_k}\right|,\hfill & \widetilde{P_k}<0\hfill \end{array}\right.$$

In these expressions, the summation includes all values of the time-series $t_j,j=1,\cdots ,M$ obtained as a result of modeling the ice growth rates $R(x,y,t_j)$ at these points in time. We similarly introduce the average ice growth rate during the periods of the positive and negative phases due to dynamic factors (advection) $\overline{A_k^{+}}$ and $\overline{A_k^{-}}$ and thermodynamic transitions (freezing) $\overline{F_k^{+}}$ and $\overline{F_k^{-}}$. For this, instead of $R(x,y,t_j)$, we will use the rates of dynamic $A(x,y,t_j)$ and thermodynamic $F(x,y,t_j)$ ice growth obtained as a result of modeling.

Figure 5 shows the fields of $\overline{R_k^{+}}$ and $\overline{R_k^{-}}$, as well as their difference $\overline{R_k^{+}}-\overline{R_k^{-}}$ for the first three components of the EOF decomposition ($k=1,2,3$), which turned out to be non-degenerate. The difference $\overline{R_k^{+}}-\overline{R_k^{-}}$ characterizes the range of changes in the ice growth rate during the transition from the negative to the positive phase.

For the convenience of comparing the contributions of different modes for several selected regions shown in Figure 5g–i,

Table 2. Standard deviations of the ice growth rate and the amplitude (half the range) of changes in the ice growth rate associated with transition from negative to positive phase of one of the EOF modes.

Regions(EOF Mode)	A(3)	B(1/2/3)	CA(1)	EL(2)	FO(2)	LK(1)
1. Standard deviations	0.97	1.31	0.89	1.19	0.69	1.07
of ice growth rate (cm/day)
including intraannual variation
2. Standard interannual deviations	0.029	0.026	0.029	0.027	0.018	0.011
of ice growth rate (cm/day)
3. Amplitude of the interanual ice	0.034	0.021/0.010/0.018	0.029	0.011	0.011	0.014
growth rate changes (cm/day)
4. Standard interannual deviations	0.029	0.091	0.027	0.049	0.024	0.041
of ice growth rate (cm/day)
due to dynamics
5. Amplitude of the interanual ice	0.025	0.071/0.097/0.094	0.027	0.025	0.014	0.008
growth rate changes (cm/day)
due to dynamics
6. Standard interannual deviations	0.022	0.092	0.018	0.043	0.013	0.041
of ice growth rate (cm/day)
due to thermodynamics
7. Amplitude of the interanual ice	0.009	0.050/0.106/0.076	0.002	0.014	0.003	0.006
growth rate changes (cm/day)
due to thermodynamics

gives the standard deviations of the ice growth rate and the amplitude (half the range) of changes in the ice growth rate associated with one of the EOF modes. According to this table’s results, the first thing to note is the smallness of interannual changes compared to seasonal ones; interannual changes make only 2–3%.

Secondly, when comparing the second and third lines of the table, we draw attention that a significant fraction of interannual variability in the selected regions is explained precisely by the variability of the principal components of the corresponding EOF modes. The only exception is the EL region, where the amplitude is less than half the standard deviation.

Third, among the selected areas, in some (A, B, EL, and LK), the fraction of variability associated with the dynamic factor differs little from the fraction of variability associated with thermodynamics (see lines 4 and 6). In contrast, in the rest (CA and FO), the dynamic factor prevails. However, irrespective of this, the EOF modes’ variability significantly better explains the dynamic interannual variability (see lines 5 and 4) rather than the thermodynamical one (see lines 7 and 6). It is evident since the EOF modes are based on the surface wind, which primarily affects the ice dynamics. An exception is a coastal part of the Beaufort Sea (B), where both the dynamic and thermodynamic fractions are equally successful, in particular, for the “dipole” and “Atlantic” modes. We should also pay particular attention to the Laptev and the Kara seas (LK). In this region, the variability of the dynamic and thermodynamic components is much greater than the total variability of the ice growth rate. This peculiarity is possible only if both components are interdependent so that the opposition of the other compensates for the strengthening of one of the factors. Another thing is that the “oceanic” mode explains only about a fifth or less of each component variability, but thoroughly explains the total variability (see lines 3 and 2). It means that the reason for such variability of the components in LK is not explicitly related to the wind effect of “oceanic” mode, but the general trend, apparently, still depends on it.

4.1. “Oceanic” Mode

A characteristic feature of the “oceanic” mode is ice growth in the central Arctic and the region of the East Siberian Sea (Figure 5g). Moreover, the mode’s negative phase is almost opposite to the positive phase, except for the Chukchi Sea and the Fram Strait. The decrease in ice prevails in both cases (Figure 5a,d). Besides, ice growth in the Beaufort Sea coastal part in the negative phase significantly exceeds its decrease in the positive phase. A similar situation is in the Laptev and the Kara seas. Because of this, in these areas, changes in the rate of ice growth during the transition from the negative to the positive phase are opposite to changes in the central region. The range of changes in the ice growth rate is most significant in the three areas highlighted in Figure 5g: area of positive values (that is, coinciding with the direction of change of the principal component) is in the central Arctic (CA), and two areas of negative values are in the coastal zone of the Beaufort Sea (B) and the area of the islands of the Laptev Sea and the Kara Sea (LK).

Figure 6 shows the time changes of the average growth rate in the CA region (Figure 6a) and, similarly, the ice decay rate in regions B (Figure 6b) and LK (Figure 6c). The correlation coefficients of these time-series with changes in the principal component are 0.549, 0.289, and 0.451. However, Figure 6d shows that the maximum correlation occurs in the presence of a time lag. For the CA region, the maximum correlation of 0.706 is achieved provided that the changes in the average ice growth rate in this region are two years ahead of the changes in the principal component, the lead time for the Beaufort Sea (B) is three years, with a correlation of 0.509, and for the Laptev and Kara Seas (LK) correlation will be equal to 0.474 when time lag is a year.

To analyze the factors contributing to the growth or decrease of ice, we first note that the dynamic factor has a dominant influence on the formation of such a response to the “oceanic” mode. Figure 7a,d,g show the dynamic component of ice growth: $\overline{A_k^{+}}$, $\overline{A_k^{-}}$ and $(\overline{A_k^{+}}-\overline{A_k^{-}})$. The same conclusion follows from Table 2 when considering amplitudes (half the ranges) of dynamic and thermodynamic components (0.071 vs. 0.050 cm/day in area B, 0.027 vs. 0.002 in CA, 0.008 vs. 0.006 in LK). The surface wind structure in the positive phase of the “oceanic” mode suggests that a positive anomaly in the central Arctic is due to the Ekman transport’s strengthening under the influence of anomalous winds. Due to this, ice forming in the marginal seas moves to the center of the anticyclonic movement. In the negative phase, on the contrary, the Ekman transport towards the central part weakens and, accordingly, the amount of ice in the coastal regions increases. As a result, anomalies are formed in the Beaufort Sea region and along the island line in the Asian part.

Thus, fluctuations of the “oceanic” mode lead to a seesaw in the ice field between the region of the central deep-water part of the Arctic and the East Siberian Sea, and the rest of the marginal seas.

4.2. “Dipole” Mode

The structure of the prevailing winds during the formation of the “dipole” mode (Figure 7h) is such that due to their action in the positive phase, there is a dynamic increase in ice in the Canadian basin from the Chukchi Sea side (Figure 7b). On the one hand, this increase is due to the Ekman transport of ice from the marginal seas of the Siberian shelf. On the other hand, it is due to the convergence of the wind field itself in this region. At the same time, there is no noticeable influence of the Ekman transport towards the islands of the Canadian Arctic Archipelago and northern Greenland. It is because the ice here is already quite cohesive. However, one can note a significant decrease in ice in the coastal part of the Beaufort Sea, which moves under the influence of wind towards the Fram Strait, where a positive anomaly is formed.

From the opposite side of the Arctic, under the Ekman transport, ice moves from the Kara and Barents Seas towards the open Arctic and concentrates at the windward part of the islands of the Severnaya Zemlya and Franz Joseph Land. In the negative phase, the effect of the wind is directly opposite (Figure 7e). The final picture for the difference in these trends (Figure 7h) is similar to the distribution in the positive phase.

Thermodynamic factors in most areas are the opposite of dynamic ones. The accumulation of ice in the open part of the Chukchi Sea in the positive phase leads to a violation of the thermodynamic equilibrium in this region and, as a result, to more intense melting of ice (Figure 8b,h). The removal of ice from the coastal part of the Beaufort Sea and the East Siberian Sea leads to the formation of young ice in these areas. The accumulation of ice on the windward side of the islands of Severnaya Zemlya and Franz Josef Land is compensated for by their melting, and the general increase in ice removal by the East Greenland current leads to their accelerated melting along the route of this current.

As a result of the action of two anti-directional factors, the following picture is formed (Figure 5b,e,h) during the transition from the negative phase to the positive:

• due to dynamic factors ice accumulates near the Arctic side of the Fram Strait (FO) (Figure 7e), while the thermodynamic factor has a predominantly neutral effect (Figure 8e);
• there is an increase in ice in the coastal part of the Beaufort Sea (B), associated with its excessive melting (Figure 8e) in the negative phase, exceeding its dynamic growth (Figure 7e);
• ice is decreasing in the East Siberian Sea and the open part of the Laptev Sea (EL), due to its dynamic entry into this region (Figure 7e) in the negative phase, not compensated by thermodynamic melting (Figure 8e).

The amplitudes of dynamic and thermodynamic components of ice growth changes in Table 2 also support the above items.

However, estimates of the correlation coefficient of the ice growth/decrease rate in these three regions do not show a stable connection with changes in the principal component of the “dipole” mode (Figure 9). For the FO region, the correlation coefficient is only 0.275, for the B region, it is 0.147, and for the EL region, it is 0.176. Nevertheless, according to Figure 9d, the maximum value of the correlation coefficient is achieved in the presence of a time lag, so that changes in the ice growth rate in the FO and EL regions anticipate changes in the principal component by two years, and in region B they are delayed by three years. Moreover, the maximum values of the correlation coefficients in the FO, B, and EL regions are then 0.518, 0.34, and 0.271, respectively. In the first case, the correlation is significant. In the second, despite the smallness, the corresponding coefficients for the dynamic and thermodynamic components are equal to −0.755 and 0.787, which indicates the simultaneous intensification of the processes of ice freezing and removal from this region in region B. In contrast, a small final value arises due to opposite directions of both processes. The final result of their actions depends on other factors. The dependence of the ice growth rate in the EL region on the variability of the “dipole” mode remains in doubt since the value of the correlation coefficient is much lower than the accepted critical value (0.5), and the expansion into dynamic and thermodynamic components gives values of the coefficients 0.262 and −0.187, respectively, and does not change the big picture.

4.3. “Atlantic” Mode

Changes in the ice growth rate due to the variability of the “Atlantic” mode are shown in Figure 5c,f,i. The greatest changes take place in the central Arctic and the sector of the East Siberian Sea (A), while changes in the opposite direction occur in the coastal region of the Beaufort Sea (B) and off the northern coast of Greenland. The main features of this variability coincide, except for details, with variability due to the “oceanic” mode. However, unlike the “oceanic” mode, a positive anomaly in region A is formed due to a significant increase in ice in the central Arctic during the positive phase (Figure 5c) and an equally significant decrease in its amount in the region of the East Siberian Sea during the negative phase (Figure 5f). In contrast, in the case of the “oceanic” mode, these regions act more synchronously (Figure 5a,d). The slowdown in growth in the Laptev and Kara seas is less pronounced than in the “oceanic” mode.

To analyze the components of the ice growth rate in these regions, we first note that the growth of ice in the Central Arctic in the positive phase and its decrease in the East Siberian Sea in the negative phase both coincide with the trends caused by the ice motion (Figure 7c,f). The decrease in ice in region B also occurs due to dynamic transport, because, judging by Figure 7c,f and Figure 8c,f, the trends in dynamic and thermodynamic growth in this region are opposite and, nevertheless. As a result, we have a general trend similar to the dynamic one (Figure 7c,f).

The average ice growth rate over region A correlates with the principal component of the “Atlantic” mode with a coefficient of 0.645 (Figure 10a). However, from Figure 10c, it follows that the greatest correlation occures with a two-year lag marked in Figure 10a with a dashed line, and the maximum value then is 0.86. The correlation coefficient for region B is only 0.272 (Figure 10b). However, from Figure 10c, it follows that the greatest correlation of 0.456 occures if these time-series are shifted by two years relative to each other. Besides, the correlations of the dynamic and thermodynamic components are more significant, and even in the absence of a shift, they are equal to 0.563 and −0.512, respectively.

Thus, the direct influence of the “Atlantic” mode turned out to be also prominent and statistically significant. Besides, it is not limited to changes in the ice growth rate in the Arctic. For example, Ref. [44] states that the negative phase (in our understanding), corresponding to the anticyclonic wind circulation in the Greenland and Norwegian seas, contributes to the atlantification of the Barents and Kara seas and an additional influx of Atlantic warm waters into the Arctic Ocean basin.

To analyze the influence of the “Atlantic” mode, we will use the heat flux through the Fram Strait and the Barents opening and also the general transport of Atlantic waters through these sections in comparison with the time variations of this mode principal component, taken with the opposite sign. As a signature of Atlantic waters at these sections, we will take, (as in [45]) temperature. We will consider waters with temperatures above 2 ${}^{\circ }$C as Atlantic waters. Following this, the transport of Atlantic waters through the vertical section L equal to

$$\mathrm{A}\mathrm{t}\mathrm{l}. \mathrm{T}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}={\int }_L{\int }_{-H}^0\left(\mathbf{u}·{\mathbf{n}}_l\right)\theta (T-2)dzdl,$$

and the corresponding heat flux

$$\mathrm{A}\mathrm{t}\mathrm{l}. \mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t} \mathrm{F}\mathrm{l}\mathrm{u}\mathrm{x}={\int }_L{\int }_{-H}^0{c}_p\rho \left(T-2\right)\left(\mathbf{u}·{\mathbf{n}}_l\right)\theta (T-2)dzdl,$$

where $dl$ is the section element of L, H is the depth of the ocean, ${\mathbf{n}}_l$ is the unit normal vector to L directed to the Arctic, $\mathbf{u}$ is the current velocity vector, $\theta$ is the Heaviside step function, T is the water temperature, $\rho$ is its density and $c_p$ is its heat capacity. The corresponding time-series are shown in Figure 11. The inflow of Atlantic waters at the entrance to the Arctic Ocean has a significant linear relationship with the third mode principal component in the presence of a three-year time lag (Figure 11c). The corresponding correlation coefficient is 0.61 (in the absence of a lag, it is only 0.269), and the standard interannual deviation is about 0.45 Sv. Moreover, the inflow through the Barents opening is more related to the third mode, and the correlation coefficient in the presence of a 5-year lag reaches 0.5. The inflow through the Fram Strait has no significant connection with the variability of the third mode.

However, the heat flux from the Atlantic depends not only on the amount of incoming water but also on their heat content. The temperature of the incoming Atlantic waters is also associated with the variability of the third mode (or NAO variability). However, compared to the flow rate, this value is more inert, and the response to changes in atmospheric circulation occurs with a large delay since advection of warmer waters requires extra time. As a result, the optimal time lag for the maximum correlation increases by about one year, and the correlation coefficient decreases to the level of 0.51 for the total heat influx from the Atlantic, remaining significant, and to 0.448 for the heat influx through the Barents Sea.

5. Discussion

As noted earlier, the maximum correlation between the ice growth rate in the most sensitive areas and the corresponding principal component occurs if the latter is shifted relative to the first in time. Moreover, a positive lag seems quite natural if we assume that changes in the ice growth rate occur due to the conditions formed by the atmosphere. However, in several cases, there is a two-three-year lead, which requires special attention and explanation. Possible explanations may include the following:

• dependence of the polar atmospheric circulation on the nature of ice distribution in the Arctic;
• dependence of the polar atmospheric circulation and the nature of the ice distribution in the Arctic on something third, and ice reacts to a change in this “third” earlier;
• some features of ice reaction to periodic wind forcing.

Undoubtedly, the first involves a broader study within the framework of a model (or series of models) of ESM class [46] with a full set of feedbacks. In our case, we obtained our results with a fixed atmospheric forcing, which although obtained in interaction with the ice dynamics, is different from the dynamics of our SibCIOM model. We can only hope that it does not differ much from our results.

We can interpret the second assumption in a more global sense, for example, the cause may be such global processes as ENSO, changes in solar activity, long-term internal oscillations of the ocean, or some other causes. This assumption, like the previous one, is beyond the scope of this particular study.

Finally, let us explain what we mean by the features of the ice reaction to periodic wind forcing. In this case, of course, we mean variations with a characteristic time scale of several years. The fact that the time-series associated with a dynamic change in the ice field is phase-shifted relative to the wind forcing indexes, i.e., maximums and minimums are reached by about two to three years earlier, can be explained by the fact that the ice reaches its equilibrium state earlier than the wind starts to change in the opposite direction. At the same time, a further increase in the wind force has no longer a significant effect on ice, since ice may already have been sufficiently grouped in certain areas and physically cannot reach any higher concentration. In models such as CICE, a further increase in wind exposure is offset by internal stresses in the ice and an increase in the effective viscosity of the ice. The mechanism of such a reaction can be understood mathematically using an idealized example of a linear oscillator with friction.

Let the ice thickness be the sum of the two terms $h=h_0+\chi$, where $h_0$ is the thickness that results from the balance of melting and freezing, about which we assume that it has reached a certain stationary state, and $\chi$ represents a deviation from this state. Let us assume that this deviation is a quantity oscillating with friction under the action of periodic perturbation. As a perturbation, let us imagine an oscillation of the principal component of some EOF decomposition mode (although an oscillation of any other origin can be assumed). For the value $\chi$ in the general case of linear systems, one can write the differential equation in the form

$$\ddot{\chi }+2\nu \dot{\chi }+{\omega }_0^2\chi =f_{0}cos\omega t,$$

where $\nu$ is the coefficient of resistance to the process occurring at a speed of $\dot{\chi }$, ${\omega }_0$ is the natural frequency of the oscillator, which characterizes the strength of the resistance to the system deviation from its equilibrium state, $\omega$ is the frequency of forcing disturbance.

The assumption that ice can be represented as a linear oscillator with friction is rather crude, and it is important to understand that the fundamental issue is to take into account nonlinear and, even to some extent, stochastic interactions. However, if we assume that there is some more or less thermodynamically stable state of ice, then after a small disturbance, the system will tend to return to this state with a characteristic time $\tau$, which will correspond to the natural frequency ${\omega }_0=2\pi /\tau$. It is easy to assume that the system will pass through the equilibrium state because of its inertia and will experience a perturbation opposite to the initial one. Therefore, it may experience some oscillations. In this case, it is also appropriate to consider some attenuation of the process with a characteristic time of $1/\nu$.

We assume that the perturbation $f_{0}cos\omega t$ is associated with the Ekman component of the wind stress. The solution to this equation can be represented as

$$\chi \left(t\right)=exp\left(-\nu t\right)\left(c_{1}cos{\omega }_{d}t+c_{2}sin{\omega }_{d}t\right)+$$

$$+\mathrm{R}\mathrm{e}\left[\frac{f_0\left({\omega }_0^2-{\omega }^2+2i\nu \omega \right)}{{\left({\omega }_0^2-{\omega }^2\right)}^2+4{\nu }^2{\omega }^2}exp\left(-i\omega t\right)\right],$$

where ${\omega }_d={\omega }_0\sqrt{1-{\nu }^2}$, and the constants $c_1$ and $c_2$ are determined from the initial conditions. If we assume that the characteristic decay time $1/\nu$ of natural oscilations is small compared with the characteristic time of external influence of $1/\omega$, i.e., $\nu \gg \omega$, then the response to such an effect can be obtained by passing to the limit at $t\to \infty$, in the end we get

$$\chi (t\to \infty )=\frac{f_0}{\sqrt{{\left({\omega }_0^2-{\omega }^2\right)}^2+4{\nu }^2{\omega }^2}}cos\left(\omega t+\varphi \right),$$

where

$$tan\left(\varphi \right)=-\frac{2\nu \omega }{{\omega }_0^2-{\omega }^2}.$$

(1)

If we believe that ice thickness is a lagged response to atmospheric forcing, then it naturally to suggest that $\varphi <0$, i.e., ${\omega }_0>\omega$. It means that time lag will be $\frac{\left|\varphi \right|}{\omega }$. In our study we considered the rate of ice growth, i.e., $\dot{\chi }$. Using previous expression we can find

$$\dot{\chi }(t\to \infty )=\frac{\omega f_0}{\sqrt{{\left({\omega }_0^2-{\omega }^2\right)}^2+4{\nu }^2{\omega }^2}}cos\left(\omega t+\varphi +\frac{\pi }{2}\right),$$

From this expression it can be seen that a phase shift similar to that detected earlier is possible if $\left(\varphi +\frac{\pi }{2}\right)>0$ or $\varphi >-\frac{\pi }{2}$. In combination with previos we finally have

$$-\frac{\pi }{2}<\varphi <0.$$

Consider the frequency spectrum of dynamic ice growth oscillations. Figure 12 shows the spectrum of this quantity averaged over the CA region (see Figure 5g). Figure 6 also shows us that for the forced oscillation of “oceanic” mode, the period is of about 15–20 years, and the optimal phase shift is of two years ahead. It means that $\left(\varphi +\frac{\pi }{2}\right)\approx 35-{45}^{\circ }$ or $\varphi \approx -(45-{55}^{\circ })$. Among the highs of Figure 12, the main one, corresponding to the annual period, could be distinguished. Besides, there are several additional highs at 6.4, 12.8, and 16.2 (in years). The latter is in a 15–20 interval estimate and most likely corresponds to the forcing frequency $\omega$. If we consider the first period (6.4 years) to be associated with the natural oscillation frequency ${\omega }_0$ then the more accurate evaluation will be $\left(\varphi +\frac{\pi }{2}\right)\approx 44.4^{\circ }$ or $\varphi \approx -(45.6^{\circ })$ and from (1) we have

$$\frac{2\nu \omega }{{\omega }_0^2-{\omega }^2}\approx 1.02.$$

Thus, we can derive an estimate for characteristic relaxation time of the ice field: $1/\nu \approx 5.9$ years. It means that relaxation time is of the same order as the natural oscillation period, and therefore this oscillation can hardly be observed.

If we consider the second period (12.8 years), which is double the first, as the natural oscillation frequency ${\omega }_0$, then another estimate for characteristic relaxation time will be $1/\nu \approx 53$ years. It means that relaxation time is much longer than the natural oscillation period, and therefore this oscillation is observable but never has been reported. From this, we conclude that the second period just corresponds to half the frequency of ice oscillation and is not a separate new oscillation.

Thus, assuming that the equation of a linear oscillator with friction can describe long-term changes in ice, we considered its response to a periodic external impact associated with long-term variability of wind circulation in the surface layer of the atmosphere. As a result, we came to a seemingly not contradictory picture, provided that the ice (its thickness) on a scale of several years has its natural oscillation frequency of about six years and approximately the same relaxation period. However, in a practical sense, at such time scales, a host of other faster and more intensive processes are observed, and in this case, even seasonal changes are quick. Therefore, it is not possible to reliably establish the validity or failure of all these estimates. In our next study, we are planning to check this approach by performing a model simulation with ideal periodic atmospheric forcing. So far, this approach can be used as a plausible hypothesis when explaining that the maximum correlation between the ice growth rate and changes in the wind field on such time scales occurs provided that one accounts for a time shift of two to three years, and changes with ice occur before wind changes.

6. Conclusions

Based on the surface wind EOF decomposition, we obtained three non-degenerate modes of the atmospheric wind forcing in the Arctic. The first (“oceanic”) mode is associated with the cyclonic or anticyclonic type of circulation in the ocean and the AOO index. The second (“dipole”) helps accelerate or slow down the Transpolar Drift. The third (“Atlantic”) weakens or intensifies the cyclonic gyre in the Northern Atlantic seas (Greenland and Norwegian).

Analyzing the rate of ice change in the Arctic associated with the “oceanic” mode, we concluded that this mode’s variability leads to the formation of a seesaw in the ice field between two regions. From the one side, it is the region of the central deep-water part of the Arctic, including the East Siberian Sea, and from the other side, it is an area including all other marginal seas. The “dipole” mode is most associated with an increase/decrease in the ice thickness at the Arctic exit near the Fram Strait and the formation of the so-called “ice factory” in the Beaufort Sea coastal regions in the positive phase of this mode. The direct influence of the “Atlantic” mode on the ice field formation is also prominent. Moreover, there is a significant relationship between the variability of this mode and the arrival of Atlantic waters with a high heat content into the Arctic through the Barents opening, which creates preconditions for reduced ice formation in this region, corresponding to the atlantification trends in the Barents Sea.