On the Annual and Semi-Annual Components of Variations in Extent of Arctic and Antarctic Sea-Ice

Abstract

In this paper, the 1978–2022 series of northern (NHSI) and southern (SHSI) hemisphere sea ice extent are submitted to singular spectral analysis (SSA). The trends are quasi-linear, decreasing for NHSI (by 58,300 km${}^2$/yr) and increasing for SHSI (by 15,400 km${}^2$/yr). The amplitude of annual variation in the Antarctic is double that in the Arctic. The semi-annual components are in quadrature. The first three oscillatory components of both NHSI and SHSI, at 1, 1/2, and 1/3 yr, account for more than 95% of the signal variance. The trends are respectively 21 (Antarctic) and 4 times (Arctic) less than the amplitudes of the annual components. We next analyze variations in pole position (PM for polar motion, coordinates $m_1$, $m_2$) and length of day (lod). Whereas the SSA of the lod is dominated by the same first three components as sea ice, the SSA of the PM contains only the 1-yr forced annual oscillation and the Chandler 1.2-yr component. The 1-yr component of NHSI is in phase with that of the lod and in phase opposition with m₁, while the reverse holds for the 1-yr component of SHSI. The semi-annual component appears in the lod and not in $m_1$. The annual and semi-annual components of NHSI and SHSI are much larger than the trends, leading us to hypothesize that a geophysical or astronomical forcing might be preferable to the generally accepted forcing factors. The lack of modulation of the largest (SHSI) forced component does suggest an alternate mechanism. In Laplace’s theory of gravitation, the torques exerted by the Moon, Sun, and planets play the leading role as the source of forcing (modulation), leading to changes in the inclination of the Earth’s rotation axis and transferring stresses to the Earth’s envelopes. Laplace assumes that all masses on and in the Earth are set in motion by astronomical forces; more than variations in eccentricity, it is variations in the inclination of the rotation axis that lead to the large annual components of melting and re-freezing of sea-ice.

1. Introduction

In a series of previous papers, we have applied singular spectrum analysis (SSA) to long time series of geophysical observables. We have analyzed series for which several decades to centuries of data were available, allowing for exploration of their periodic or quasi-periodic components on time scales from days to a full century. We have successively analyzed polar motion and length of day [1,2,3], the oceanic climate indices MJO, PDO, ENSO, etc. [4,5,6], and global surface temperatures [7,8], and contributed to the understanding of the solar–terrestrial relationship of sunspot numbers via ISSN [9,10]. The determination of a comprehensive range of pseudo-periodic components of ISSN has even allowed us to propose a prediction of the ongoing solar cycle, Cycle 25. In the present paper, we extend the same analysis to an important component of the climate system, namely, variations in the surface extent of sea ice in the Arctic and Antarctic.

In previous papers, on which the present study is partly based [2,3], we recalled the importance of Laplace’s treatise on celestial mechanics [11], in particular the system of equations, later named after Euler and Liouville, that express the conservation of kinetic momentum [12] and govern the declination and inclination of the rotation axis of any planet in revolution about the Sun (see Appendix A and Appendix B). Laplace’s paradigm can be summarized in two statements: (1) all masses on the surface or internal to a planet are set in motion by astronomical forces (particularly from the Sun, Moon, and Jovian planets) and (2) it is sufficient to know the polar motion or its time derivative (the lod) in order to have access to the periods of all moving masses on the Earth’s surface or inside of it. This is why the revolution periods of planets (or their commensurate combinations) and luni-solar tides have been found in many observed time series of geophysical phenomena, ranging from hours to millions of years [1,2,5,7,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30].

Most recently, Courtillot et al. [10] have shown that the sunspot number ISSN, which is a proxy for solar activity, can be decomposed as a series of oscillations with periods similar to or commensurate with those of the Jovian planets [14]. They were able to reconstruct more than 90% of the variance of the original time series [10]. Similar results have been obtained for polar motion, with 90% of the variance being carried by the ephemerids of the Jovian planets [2]. Pioneering work by Milankovic [13] and more recently by Laskar et al. [17] explains how, over very long time scales, these planets influence the obliquity, precession, and eccentricity of our planet’s orbit and rotation axis, and therefore, among other consequences, its climate.

Planets exert gravitational forces (attraction) on all celestial bodies; because of its mass, the Sun is the leading attractor, exerting a force F of 3.57 × 10${}^{22}$ kg.m.s${}^{-2}$ on the Earth. However, because it is close to the solar system’s center of gravity, its contribution to the budget of angular momenta (F*Sun–Earth distance*Earth’s revolution period = 1.68 × 10${}^{41}$ kg.m${}^2$.s${}^{-1}$) is one to two orders of magnitude smaller than that of the Jovian planets. In a heliocentric reference frame, Jupiter’s angular momentum is 1.93 × 10${}^{43}$ kg.m${}^2$.s${}^{-1}$, Saturn’s is 7.82 × 10${}^{42}$ kg.m${}^2$.s${}^{-1}$, Neptune’s is 2.50 × 10${}^{42}$ kg.m${}^2$.s${}^{-1}$, and that of Uranus is 1.69 × 10${}^{42}$ kg.m${}^2$.s${}^{-1}$.

We present the data for the extent of sea ice in the Northern and Southern hemispheres (NHSI and SHSI, respectively) in Section 2, then their spectral components, using both Fourier analysis and SSA, in Section 3. We provide complementary results on polar motion and lod in Section 4, then compare the annual components of polar motion and lod with those of the sea ice extent in Section 5. Finally, we discuss a hypothesis based on these results and close with a conclusion in Section 6.

2. Available Data on the Extent of Sea Ice in Both Hemispheres

Very precise daily series of sea ice extent in the Arctic (NHSI) and Antarctic (SHSI) are available starting in 1978 ([31,32], https://nsidc.org/data/seaice_index/data-and-image-archive, accessed on 2 April 2022). The dataset from which the NHSI and SHSI sea ice areas are calculated consists of sea ice concentration maps derived from the radiance obtained from microwave radiometers on a suite of satellites. These observations began with the Nimbus 7 Scanning Multichannel Microwave Radiometer (SMMR), which operated from 1978 through 1987, followed by the Defense Meteorological Satellite Program (DMSP) series of F8, F11, F13, and F15 Special Sensor Microwave Imagers (SSMIS) and the F17 Special Sensor Microwave Imager Sounder (SSMIS). The F8 satellite operated from 1987 through 1991, the F11 from 1991 through 1995, the F13 from 1995 through 2007, and the F17 from 2008 through 2010 [31,33,34,35]. After 2010 and up to 2018, we refer to [36], quoting [34,35].

The dataset is generated using the Advanced Microwave Scanning Radiometer—Earth Observing System (AMSR-E) Bootstrap Algorithm with tie-points varying daily. Daily (every other day prior to July 1987) and monthly data are available for both the north and south polar regions. Data are gridded on the SSM/I polar stereographic grid (25 × 25 km) and provided in two-byte integer format. The sampling from 26 October 1978 to 2 December 1987 is not as continuous as afterwards; there is only one datum every other day from 26 October 1978 to 20 August1987, one per day from 20 August 1987 to 2 December 1987, then a gap from 2 December 1987 and 13 January 1988, back to one point per day from 13 January 1988 to 19 February 2022, a two day gap from 19 February 2022 to 22 February 2022, and finally full sampling (one per day) until 21 March 2022. We worked with the dataset as described above, starting in 26 October 1978 and ending in 21 March 2022 (see https://nsidc.org/data/G02135/versions/3 for more details, accessed on 2 April 2022). The data are shown in Figure 1a, with NHSI in blue and SHSI in red.

A number of features are immediately apparent from this figure: annual variations are dominant, their amplitude in the Antarctic is double that in the Arctic, and they are in phase opposition and modulated, with a flat trend for SHSI, a decreasing trend for NHSI, and a sub-decadal modulation of the envelopes of both. Figure 1b shows the Fourier spectra of the two series with successive enlargements; it can be seen that a fundamental annual period is accompanied by seven harmonics ranging from 1/2 to 1/8 of a year.

3. SSA Analysis of NHSI and SHSI Data on Sea Ice Extent

As in a number of previous papers (e.g., [2,5,10]), we have submitted time series to iterative Singular Spectrum Analysis (iSSA), and now do the same for the sea ice extent. We refer the reader to these papers and to the Golyandina and Zhigljavsky’s book [37] for the SSA method, to [38] for the properties of the Hankel and Toeplitz matrices that it uses, and to [39] for the singular value decomposition algorithm SVD).

Let us begin by considering a discrete (non-zero) time series (${\mathcal{X}}_N$) of length N (N > 2): ${\mathcal{X}}_N=(x_1,\dots ,x_N).$

Step 1: Embedding Step

${\mathcal{X}}_N$ is divided into K segments of length segments of length L in order to build a matrix $\mathbf{X}$ with dimensions $K\times N$, where $K=N-L+1$ conditions our decomposition. This is the first “tuning knob”. Integrating $\mathbf{X}$ yields a Hankel matrix

$$\mathbf{X}=\left(\begin{array}{cccc}{x}_1& x_2& x_3\cdots & x_K\\ x_2& x_3& x_4\cdots & x_{K+1}\\ x_3& x_4& x_5\cdots & x_{K+2}\\ ⋮& ⋮& ⋮\ddots & ⋮\\ x_L& x_{L+1}& x_{L+2}\cdots & x_N\end{array}\right)$$

(1)

Embedding, which is the first step in an SSA, consists of projecting a one-dimensional time series in a multidimensional space of series ${\mathcal{X}}_N$ such that vectors $X_i={(x_i,\dots ,x_{i+L-1})}^t$ belong to ${\mathcal{R}}^L$, where $K=N-L+1$. The parameter that controls the embedding is L, the size of the window used for analysis, which is an integer between 2 and $N-1$. The Hankel matrix has a number of properties with regard to symmetry; the dimension of its transpose ${\mathbf{X}}^t$, called the trajectory matrix, is K. Embedding is a compulsory step in the analysis of nonlinear series. Formally, it consists of the empirical evaluation of all pairs of distances between two offset vectors, which are delayed (lagged) in order to calculate the correlation dimension of the series.

Step 2: Decomposition in Singular Values (SVD)

SVD of a non-zero trajectory matrix X with dimensions $L\times K$ takes the shape

$$\mathbf{X}=\sum _{i=1}^d\sqrt{{\lambda }_i}{U}_i{V}_i^t$$

(2)

where the eigenvalues ${\lambda }_i$ ($i=1,\dots ,L$) of matrix $\mathbf{S}={\mathbf{XX}}^T$ are arranged in order of decreasing amplitude. Eigenvectors $U_i$ and $V_i$ are provided by

$$V_i={\mathbf{X}}^T{U}_i/\sqrt{{\lambda }_i}$$

(3)

Here, the $V_i$ form an orthonormal basis and are arranged in the same order as the ${\lambda }_i$. Let ${\mathbf{X}}_i$ be a part of matrix $\mathbf{X}$ such that

$${\mathbf{X}}_i=\sqrt{{\lambda }_i}{U}_i{V}_i^t.$$

(4)

Embedding matrix X can then be represented as a simple linear sum of elementary matrices X${}_i$. If all eigenvalues are equal to 1, then decomposition of X into a sum of unitary matrices is as follows:

$$\mathbf{X}={\mathbf{X}}_1+{\mathbf{X}}_2+\dots +{\mathbf{X}}_d$$

(5)

As d is the rank of X ($d=\mathrm{rank}\mathbf{X}=max\left\{i\right|{\lambda }_i>0\}$), SVD allows us to write X as a sum of d unitary matrices defined in a univocal way.

We now discuss the nature and characteristics of the embedding matrix. Its rows and columns are subseries of the original time series (or signal). The eigenvectors $U_i$ and $V_i$ have a time structure, and can be considered as a representation of temporal data. Let X be a suite of L lagged parts (X and ${\mathbf{X}}_1$,…, ${\mathbf{X}}_K$) of the linear basis of its eigenvectors. If we let

$$Z_i=\sum _{i=1}^d\sqrt{{\lambda }_i}{V}_i$$

(6)

with $i=1,\dots ,d$, then relation (4) can be written as

$$\mathbf{X}=\sum _{i=1}^d{U}_i{Z}_i^t$$

(7)

that is, for the jth elementary matrix,

$$X_j=\sum _{i=1}^d{z}_{ji}{U}_i$$

(8)

where $z_{ji}$ is a component of vector Z${}_i$. This means that vector Z${}_i$ is composed of the ith components of vector X${}_j$. In the same way, if we let

$$Y_i=\sum _{i=1}^d\sqrt{{\lambda }_i}{U}_i$$

(9)

then for the transposed trajectory matrix we obtain

$$X_j^t=\sum _{i=1}^d{U}_i{Y}_i^t$$

(10)

which corresponds to a representation of the lagged vectors in the orthogonal basis ($V_1$,…,$V_d$). It is apparent why SVD is a very good choice for the analysis of the embedding matrix, because it provides two different geometrical descriptions.

Step 3: Reconstruction

As we have seen, ${\mathbf{X}}_i$ matrices are unit matrices, and we can “re-group” these matrices into a physically homogeneous quantity. This is the second “tuning knob” of SSA. In order to regroup the unit matrices, we first divide the set of indices {i${}_1$,…,i${}_d$} into m disjoint subsets of indices $\{I_1,\dots ,I_m\}$.

Let I be the grouping of p indices of $I=\{i_1,i_2,\dots ,i_p\}$; because (5) is linear, the resulting matrix X${}_I$ that regroups indices I can then be written as

$${\mathbf{X}}_I={\mathbf{X}}_{I1}+{\mathbf{X}}_{I2}+\dots +{\mathbf{X}}_{Im}$$

(11)

This step is called regrouping the eigentriplets ($\lambda$, U and V). In the limit case $m=d$, (11) becomes exactly (5), and we again find the unit matrices.

Next, we need to associate pairs of eigentriplets. This means separating the additive components of a time series. First, we must consider the concept of separability.

Let $\chi$ be the sum of two time series ${\chi }^{\left(1\right)}$ and ${\chi }^{\left(2\right)}$ such that $x_i=x_i^{\left(1\right)}+x_i^{\left(2\right)}$ for any $i\in [1,N]$. Let L be the analyzing window (with fixed length), and let X, X${}^{\left(1\right)}$, and X${}^{\left(2\right)}$ be the respective embedding matrices of series $\chi$, ${\chi }^{\left(1\right)}$, and ${\chi }^{\left(2\right)}$. These two subseries are separable (even weakly) in Equation (5) if there is a collection of indices $\mathcal{I}\subset \{1,\dots ,d\}$ such that ${\mathbf{X}}^{\left(1\right)}={\sum }_{i\in \mathcal{I}}{\mathbf{X}}_i$, respectively, if there is a collection of indices such that ${\mathbf{X}}^{\left(1\right)}={\sum }_{i\notin \mathcal{I}}{\mathbf{X}}_i$.

In a case in which separability exists, the contribution of ${\mathbf{X}}^{\left(1\right)}$ (for instance) corresponds to the ratio of the associated eigenvalues (${\sum }_{i\in \mathcal{I}}{\lambda }_i$) to the total eigenvalues (${\sum }_{i=1}^d{\lambda }_i$).

Thus, regrouping SVD components can be summarized by their decomposition into several elementary matrices, the structure of which must be as close as possible to a Hankel matrix of the initial trajectory matrix (this is true on paper; however, in reality things are much more difficult).

Step 4: The Diagonal Mean (or Hankelization) Step

The final step involved returning to the data space in order to calculate time series with length N associated with submatrices ${\mathbf{X}}_I$. Let Y be a matrix with dimension $L\ast K$; then, for each element $y_{ij}$ we have $1⩽i⩽L$ and $1⩽j⩽K$. Let $L^{*}$ be the minimum and $K^{*}$ be the maximum; in addition, one always has $N=L+K-1$. Finally, let $y_{ij}^{*}=y_{ij}\mathrm{if}L<K;\mathrm{otherwise},\mathrm{let}{y}_{ij}^{*}=y_{ji}$. The diagonal average applied to the kth index of time series y associated with matrix Y provides us with

$$y_k=\left\{\begin{array}{ccc}\frac{1}{k}& {\displaystyle \sum _{m=1}^k}{y}_{m,k-m+1}^{*}& 1⩽k⩽L^{*}\\ \frac{1}{L^{*}}& {\displaystyle \sum _{m=1}^{L^{*}}}{y}_{m,k-m+1}^{*}& L^{*}⩽k⩽K^{*}\\ \frac{1}{N-K+1}& {\displaystyle \sum _{m=k-K*+1}^{N-K*+1}}{y}_{m,k-m+1}^{*}& K^{*}⩽k⩽N^{*}\end{array}\right.$$

(12)

Relation (12) corresponds to the mean of the element on the anti-diagonal $i+j=k+1$ of the matrix, while for $k=1$, $y_1=y_{1,1}$, for $k=2$,$y_2=(y_{1,2}+y_{2,1})/2$, etc.

As mentioned above, Step 3 is the most difficult part. Here, we have chosen the iterative SSA (iSSA) approach from among many others. Because relation (5) is linear, we can iterate the decomposition. We start with a small value of L (looking for the longest period), which we then increase until we obtain a quasi-Hankel matrix (Steps 1 and 2). We then extract the corresponding lowest-frequency component that was subtracted from the original signal. We again increase the value of L to find the next component (shortest period). The algorithm stops when no pseudo-cycle can be detected or extracted. In this way, we scan the series from low frequencies to high.

We next apply SSA in order to determine the waveforms (amplitude and phase) of the periodic or quasi-periodic components of the sea ice extent series. As seen in Figure 1, the series are dominated by the modulation of a quasi-monochromatic cycle.

For NHSI, the first three components allow for recovery of 98.7% of the total signal variance. They occur at periods of 1.000 ± 0.014 yr (95.6% of the variance), 0.499 ± 0.003 yr (2.69%), and 0.333 ± 0.001 yr (0.43%). For SHSI, the first three components permit recovery of 97.9% of the total signal variance, occurring at periods of 0.999 ± 0.014 yr (94.4%), 0.500 ± 0.004 yr (3.25%), and 0.333 ± 0.001 yr (0.25%). The remaining components have negligible contributions. In the following, we focus on the annual and semi-annual components, which are respectively displayed in Figure 2a,b.

The annual Antarctic sea ice component is only slightly modulated on a multi-decadal time scale, whereas the annual Arctic component is modulated, with significant growth. The two are in phase opposition (180 ± 1 day), as expected. For the semi-annual components, the phase difference is 34 ± 2 days, which is close to quadrature (90/2 days).

The trends are quasi-linear (not shown in the figures), decreasing for NHSI (by 58.300 km${}^2$/yr) and increasing for SHSI (by 15.400 km${}^2$/yr). The total surface change since 1978 is a gain of 0.57 × 10${}^6$ km${}^2$ for the Antarctic and a loss of 2.15 × 10${}^6$ km${}^2$ for the Arctic. These values are in good agreement with [36]; note that they are respectively 21 (Antarctic) and 4 (Arctic) times less than the amplitudes of the annual SSA components.

4. Complementary Results for Polar Motion and Length of Day

The orientation of the Earth’s rotation axis is a key variable on many time scales that modulates the insolation of any location on the Earth’s surface. It is therefore interesting to compare the spectral components of sea ice extent with the variations in polar motion (PM) and length of day (lod). These data are available from the website of the International Celestial Reference System (IERS, https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html, accessed on 2 April 2022), from which we have selected File EOP14C04. Several authors have studied the spectral content of the lod. In Le Mouël et al. [1], we used SSA to explore the range from 9.13 days to 18.6 years. Longer periods (at 18.6 yrs and 11 yrs) as well as certain harmonics, along with shorter periods at 27 days and associated harmonics, are linked to the Sun and Moon, and are discussed elsewhere [1]. In this paper, we focus on the range below 2 years (Figure 3, Fourier spectrum). Distinct peaks can be seen at 1 and 1/2 yr, as well as at 1/3 and 1/4 yr. The leading two terms amount to 13.25% and 11.06% of the total lod variance, respectively, for a total of 24.3%.

We now turn to the time series (1978–2022) of coordinates of the pole of rotation ($m_1$, $m_2$; see [2]). The Fourier analysis of the modulus m (hereafter PM, for polar motion) is shown in Figure 4. The spectrum features only two sharp peaks at 1 yr and 1.2 yr. Components with decadal to ∼80 yr (Gleissberg cycle) periods are present, though out of the picture, as is the trend known as the Markowitz drift [40,41,42]. The 1/3 yr and higher harmonics seen in NHSI, SHSI, and lod are not found in PM The ∼1.2 yr component is actually a doublet at 1.19 ± 0.00(4) and 1.20 ± 0.00(4) yr, known as the Chandler wobble [43,44]. These three components, Markowitz, Chandler and annual, respectively amount to 7.5, 40.4, and 19.8% of the signal variance, for a total of 67.7% [2].

Figure 5a,b compares the annual and semi-annual SSA components of PM and lod. The annual components are slightly modulated (<10%) and are close to pure sinusoids. As a result, they are almost perfectly correlated, with a phase difference of 41 ± 3 days between $m_2$ and lod (i.e., quadrature; Figure 5).

The prominence of an “astronomical” annual component in these time series, of course, leads one to relate it to the revolution period of the Earth around the Sun. Here, we compare the annual variations of the Sun–Earth distance D${}_{SE}$ (in au) with m₁ (Figure 6, bottom) and lod (Figure 6, top). The phase lags are respectively 65.6 ± 2.1 days for D${}_{SE}$ and $m_1$ and 152.5 ± 4.1 days for D${}_{SE}$ and lod.

5. Comparison of Annual and Semi-Annual SSA Components of Polar Motion and Length of Day vs. Variations in Sea Ice Extent

We now compare the annual components of sea ice extent NHSI (Figure 7) and SHSI (Figure 8) with those of $m_1$ and lod, in particular, their phases (

Table 1. Phase differences of the annual lines of the pairs, shown as column and line headings (in days, which in this case is very close to 1 day = 1 degree).

	D${}_{SE}$	NHSI	SHSI
lod	152.5 ± 4.1	33.0 ± 1.1	153.7 ± 7.1
$m_1$	65.6 ± 2.1	153.8 ± 4.2	11.5 ± 2.4

). In the case of the northern hemisphere (Arctic), variations in sea ice extent are almost in phase with the lod and in phase opposition with $m_1$. The phase difference is 153.8 ± 4.2 days for m${}_1$ and 33.0 ± 1.1 days for lod (see Figure 7, middle enlargement). NHSI is more strongly modulated (by 20% over 40 years, corresponding mainly to the trend) than lod. In the case of the southern hemisphere (Antarctic), variations in sea ice extent are almost in phase with m₁, and are in phase opposition with lod.

The semi-annual component of lod is compared with that of NHSI in Figure 7 (middle) and with that of SHSI in Figure 9 (top). The respective phase differences are 51.8 ± 1.9 days (N) and 19.2 ± 1.5 days (S). Recall that the semi-annual component is present in lod and is not present in $m_1$. Thanks to the Euler–Liouville system ([11]), we can integrate the lod to obtain the theoretical polar motion that should accompany it, which is shown in Figure 9 (bottom). The agreement in terms of the phase is excellent.

As shown in the previous section, these annual and semi-annual components of NHSI and SHSI are much larger (up to 20 times) than the trends observed since 1978. The mechanisms generally suggested as a forcing factor of the trends (global warming, winds; see [36,45,46,47]) must be dwarfed by a first order geophysical or astronomical forcing. Of course, the annual oscillating components of NHSI and SHSI are in opposite phases in both hemispheres (Figure 2a,b). SHSI is not significantly modulated, whereas NHSI is, which is reminiscent of the classical behavior of a forced oscillator. Polar motion indeed possesses a forced annual component, generally attributed to climate and oceanic/atmospheric sources (see, e.g., [48] or [49] chapter 7).

The lack of modulation of the largest forced component (SHSI) suggests an alternate mechanism. If we follow Laplace’s paradigm, the torques exerted by the Moon, Sun, and planets play the leading role (if not the whole) as the source of forcing (modulation) of many (if not most) geophysical phenomena (generalized Earth tides). These forces (and torques) lead to changes in the inclination of the Earth’s rotation axis, transferring stresses to the Earth’s solid and fluid envelopes, setting Earth masses, and resulting in thermal dissipation. This description is essentially that of [13] with respect to paleo-climates and ice ages; more than variations in eccentricity, it is variations in inclination of the rotation axis that lead to the large annual component of the melting and re-freezing of sea ice.

We have seen in Figure 6 the excellent correlations of lod and $m_1$ with the Sun–Earth distance D${}_{SE}$ in phase opposition and phase, respectively. It is noteworthy that the phase differences of the annual components of NHSI with $m_1$, SHSI with lod, and D${}_{SE}$ with lod are all equal to 153 days, which happens to be half of the Euler period of 306 days (Table 1).

Let us follow Lambeck’s formulation ([49], chapter 3, system 3.2.9) of the Liouville–Euler equations (see [2]):

$$\begin{array}{cc}\hfill {\displaystyle \frac{i}{{\sigma }_r}}\left({\displaystyle \frac{d\mathbf{m}}{dt}}\right)+\mathbf{m}& =\mathbf{f}\hfill \\ \hfill {\displaystyle \frac{d{m}_3}{dt}}& =f_3\left(lod\right)\hfill \end{array}$$

(13)

where $i=\sqrt{-1}$; $\mathbf{m}=m_1+i\ast m_2$ (polar coordinate). This is a linear first-order differential system. The annual oscillation of the forcing function $\mathbf{f}=f_1+i\ast f_2$ that follows sea ice extent variations should behave as $m_1$ or $m_2$, without phase changes. This is indeed what is seen in Figure 7, Figure 8 and Figure 9, where the phases remain constant. The fact that the phase differences of the polar motion and lod can be recovered from the leading annual components of sea ice strengthens previous findings (e.g., [2,3,10]), with lod and $m_1$, $m_2$ being strongly linked.

6. Discussion and Conclusions

The differences between the frequency ranges of the length of day (lod) and polar motion (PM) can be understood in the frame of Laplace’s ([11]) formulation of the Liouville–Euler system, which consists of an equation for polar axis inclination $\theta$ and another for the derivative of declination ${\displaystyle \frac{d\psi }{dt}}$, both depending only on $\theta$. The second equation, which involves a derivative operator, amounts to a high-pass filter (the derivative operator is the reason that lod contains richer higher frequency components than $m_2$), whereas the first behaves as an integrator, i.e., passes lower frequencies.

We now explore further consequences of Laplace’s formulation. We have seen that there is no semi-annual component in the polar motion (e.g., [2,3]). According to the theory of forced oscillators (e.g., [50]), the forcing frequency, (here the annual oscillation) is at the boundary between the two frequency domains. Figure 7 (center) shows that the semi-annual components of NHSI and lod are almost in quadrature (51.8 ± 1.9 days, close to 90/2 = 45 days for quadrature). Because there is a derivative operator between lod and PM, and because from Figure 5, Figure 6 and Figure 7 the forced annual oscillations may be assumed to follow variations in the Earth–Sun distance D${}_{SE}$ in lod, in PM, and finally in NHSI and SHSI with remarkably constant phase differences, any semi-annual oscillation of the lod should theoretically be present in the PM as well. As this would perturb polar motion, it must be eliminated (dissipated), which is why we find it in NHSI and SHSI and in quadrature with lod (Figure 7, bottom).

Integrating the semi-annual component of lod should provide access to the same semi-annual component of $m_1$ that does not exist, but which should be in phase with the semi-annual component of NHSI considering the behavior of the annual component of the sea ice series. The result is shown in Figure 7 (bottom) under the label “calculated $m_1$”. The phases and modulations of the two curves are the same (additionally yielding a prediction of the semi-annual component of NHSI). This is why all forced oscillations and their harmonics are found in the variations in sea ice, with very small amplitude modulations and remarkably stable phase differences between both hemispheres and with lod and PM.

The generally accepted origin of the annual oscillation in PM is the reorganization of fluid masses at the Earth surface. We may quote somewhat extensively Lambeck’s introduction ([49], chapter 7):

“The principal seasonal oscillation in the wobble is the annual term which has generally been attributed to a geographical redistribution of mass associated with meteorological causes. Jeffreys, in 1916, first attempted a detailed quantitative evaluation of this excitation function by considering the contributions from atmospheric and oceanic motion, of precipitation, of vegetation and of polar ice. Jeffreys concluded that these factors explain the observed annual polar motion, a conclusion still valid today, although the quantitative comparisons between the observed and computed annual components of the pole path are still not satisfactory. These discrepancies may be a consequence of (i) inadequate data for evaluating the known excitations functions, (ii) the neglect of additional excitation functions, (iii) systematic errors in the astronomical data, or (iv) year-to-year variability in the annual excitation functions”.

In other words, [49] explains that the annual oscillation of PM is in part caused by variations in sea ice and in part by variations in the atmospheric circulation. Regardless of whether one follows [11] or [48], the observations show remarkably constant phase differences; yet, there is no reason for the periodic components of all climate indices to have the same relative or absolute phases (e.g., [5], Figure 10a,b).

Figure 10a,b shows the Fourier spectrum of the climatic indices associated with the Arctic (AO, [51]) and Antarctic (AAO, [51]) sea ice thickness. It is common to find suggestions in the literature that the evolution of sea ice extent should be captured by its climate indices.

The AO is characteristic of exchanges in atmospheric mass between the Arctic Ocean and the zonal ring centered on 45° N. For instance, [52] have shown how changes in surface winds associated with the trends and fluctuations of AO influence sea ice motion in the center of the Arctic, and how these changes in ice motion in turn influence the thickness and concentration of sea ice and the distribution of surface temperature above the Arctic ([53,54,55]). The AAO monitors the periodic reinforcement and weakening of the circumpolar vortex, and is sometimes called the austral annular mode AAM. It is a stationary mode, unlike the Antarctic circumpolar waves ([56]) that explains a large part of Antarctic climate variability on pluri-annual scales. These two indices, AO and AAO, are widely used to explain the evolution of sea ice since 1978; yet, their spectra (Figure 10a,b) are devoid of any significant annual and semi-annual periodicities, while the NHSI and SHSI series display annual and semi-annual periodicities that amount to more than 98% of the total signal variance. Nor are annual and semi-annual periodicities found in either AO or AAO. Overall, neither annual nor semi-annual periodicities are found in either ten MJO indices or several other index series, such as PDO and ENSO ([5]). On the other hand, polar motion and length of day are global phenomena that could explain the regularity of the phase and amplitude modulations of oscillating sea ice components at 1 and 0.5 year.

Thus, Laplace’s formulation by and modern analysis of the periodic and semi-periodic components of sea ice variations are not that far from one another, although the direction of causality may have to be changed. If we follow Laplace’s view, variations in sea ice extent are a consequence of polar motion, as can be read from the Liouville–Euler equations. On the other hand, in the modern interpretation, they are consequences of wind speed variations that should be reflected by indices such as AAO and AO. We have already discussed the impact of the AAO index on the spatial variations of the annual and semi-annual components of the SLP ([57]).

The results obtained in this paper on the phases and amplitudes of the SSA components of NHSI and SHSI lead us to prefer Laplace’s formulation. This formulation allows us to understand the presence of the harmonic series (1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, and 1/8 yr) in Figure 1b. In a schematic way, variations in sea ice extent are the infinitesimal (incremental) expression of oceanic motion, a topic that requires use of the theory of fluid mechanics of turbulent flow. There is no solution of this problem in the spherical case [58,59,60,61,62,63,64], although there is one in the case of flow between two coaxial cylinders [65,66,67,68]. The perturbation in flow velocity $v_1$ can be written as

$$v_1(r,\phi ,z)=f\left(r\right)\ast e^{i(n\phi +kz-\omega t)}$$

(14)

where integer n is the order of symmetry (n=0 for axial symmetry) and k is the wavelength of the instability. Acceptable values of the annual oscillation are obtained by solving Equation (14) under the appropriate boundary conditions ($v_1$ = 0 for r = $R_1$ and r = $R_2$). The Reynolds number $Re$ can be taken as ${\displaystyle \frac{{\Omega }_1{R}_1^2}{\nu }}$ (or ${\displaystyle \frac{{\Omega }_2{R}_2^2}{\nu }}$, as ${\displaystyle \frac{R_1}{R_2}}$ and ${\displaystyle \frac{{\Omega }_1}{{\Omega }_2}}$ are given. For fixed n and k, the values of $\omega$ form a discrete suite ${\omega }_n\left(k\right)$. By analogy with Earth, the solid surface and free surface of the atmosphere are such that $v_1=0$, and the only solutions for abs($\omega$) are 1 day and 1 year. Therefore, we should try to find harmonics of these two values. A detailed study of the solution of Equation (2) in the case of co-axial cylinders is provided by [67]. The $v_1$ flows are stationary toric vortices, known as Taylor vortices, regularly placed along the cylinder generatrices. In the case of two cylinders rotating in the same sense, we find two vortices rotating in opposite directions, $2\pi /k_{cr}$, $k_{cr}$ being the critical wavelength above which $v_1$ is not 0. Thus, we know how to describe the figure of the turbulent flow between two co-axial cylinders.

In light of the above, it is no surprise to find the series of harmonics ${\omega }_n\left(k\right)$ of the annual forcing. Furthermore, the term $n\phi$ in Equation (14) forces the spatial symmetries of the flow. This is already known regarding the variations of the sea ice dipole in Antarctica (cf. [69,70,71]) and variations in the trends of global atmospheric pressure (cf. [6]). These spatial and temporal variations seem to be better explained by Taylor–Couette type flow forced by variations in polar motion (rotation), in essentially the same way as envisioned by Laplace [11].