Modeling the annual cycle of daily Antarctic sea ice extent

The total Antarctic sea ice extent (SIE) experiences a distinct annual cycle, peaking in September and reaching its minimum in February. In this paper we propose a mathematical and statistical decomposition of this temporal variation in SIE. Each component is interpretable and, when combined, gives a complete picture of the variation in the sea ice. We consider timescales varying from the instantaneous and not previously defined to the multi-decadal curvilinear trend, the longest. Because our representation is daily, these timescales of variability give precise information about the timing and rates of advance and retreat of the ice and may be used to diagnose physical contributors to variability in the sea ice. We define a number of annual cycles each capturing different components of variation, especially the yearly amplitude and phase that are major contributors to SIE variation. Using daily sea ice concentration data, we show that our proposed invariant annual cycle explains 29 % more of the variation in daily SIE than the traditional method. The proposed annual cycle that incorporates amplitude and phase variation explains 77 % more variation than the traditional method. The variation in phase explains more of the variability in SIE than the amplitude. Using our methodology, we show that the anomalous decay of sea ice in 2016 was associated largely with a change of phase rather than amplitude. We show that the long term trend in Antarctic sea ice extent is strongly curvilinear and the reported positive linear trend is small and dependent strongly on a positive trend that began around 2011 and continued until 2016.

1 Introduction

Much of the research on Antarctic sea ice variability focuses on the monthly, seasonal and interannual timescales (Parkinson and Cavalieri, 2012; Simpkins et al., 2012; Holland, 2014; Turner et al., 2015b; Hobbs et al., 2015; Holland et al., 2017). This is useful and necessary, especially if links to the larger-scale (and remote) atmospheric and oceanic forcings are to be made. However, significant aspects of the timing of the ice cycle, for example when ice advance or ice retreat begins, occur at sub-monthly scales (Stammerjohn et al., 2008; Stuecker et al., 2017; Turner et al., 2017; Schlosser et al., 2018; Meehl et al., 2019). Using daily data facilitates analysis of the daily variation in sea ice and is the springboard of this research.

The dominant or primary characteristic of Antarctic sea ice variability is its annual cycle. Satellite-observed total Antarctic sea ice extent (SIE) experiences a distinct annual cycle, peaking in September (19 million km²) and reaching its minimum in February (3 million km²) on average. In Julian days, the median minimum day is 50 and the median maximum day is 255. The growth from minimum (trough) to maximum (peak) is slower than the retreat from maximum to minimum. This is arguably the strongest seasonal cycle on the planet. The characteristics of the annual cycle that are of major interest are its amplitude and its phase. The amplitude is considered to be the difference between SIE at maximum and SIE at minimum. The phase is the timing of advance and retreat of the ice with respect to the typical annual cycle. In recent years the sensitivity of the amplitude and phase to climate change has been the subject of much study (e.g., Stammerjohn et al., 2008; Turner et al., 2017, 2015a; Parkinson, 2019).

The daily, annual cycle of SIE is traditionally calculated by simply taking the average (or the median value) for each day of the year. However, satellite-observed SIE can vary widely from day to day. Some of this variation is due to the ice growth, melting, and divergence of the ice at the ice edge and land spillover (coastal effect of mixed land/water grid cells), while some is due, for example, to transient effects of cloud and melt on the ice surface (e.g., Comiso and Steffen, 2001). A simple daily average or median includes all of these sources of variability, perhaps leading to overestimation or underestimation of the SIE. Therefore, a standard deviation (or a percentile) is often included to give some idea of the variability of the individual days around the mean for that day. While simple and transparent, this method of calculating the annual cycle produces a value that is subject to substantial variation since it is based on as few as 40 numbers (the length of the satellite-observed data time series), one for each year of recorded data, and does not include the effect of the day preceding or the day following the averaged day. It is also influenced by the pattern of missing values. Finally, it also disguises the fact that the daily annual cycle might be slowly changing phase and that the amplitude and shape of the daily annual cycle of SIE might vary. This can make it difficult to make statistically sound conclusions about variability in the data.

Our overarching aim in this research is not only to redefine the annual cycle but also to make a meaningful decomposition of the variation in the annual cycle of Antarctic SIE. We do so on the time dimension in such a way that each component can be interpreted individually, and when taken together all of the components give a complete picture of the variation in the sea ice. We consider the variation from the shortest timescale (instantaneous variation) and increasing the timescale sequentially we move through the day-to-day variation, the year-to-year (interannual) variation, and finally the longest timescale, the curvilinear trends of the multi-decadal variation. In the process, we make a number of technical contributions, most importantly to define complementary types of annual cycles that are meaningful in terms of this decomposition and also to the representation of volatility. We have deliberately chosen (time) dimensions based on their interpretability rather than solely statistical efficiency concerns. For example, the amplitude and phase components of the decomposition are much more interpretable than simple spectral components.

We begin by presenting a stochastic model for the sea ice extent that allows the annual cycle to be defined in flexible ways. This model can represent the real variability in SIE and reduces the contribution from the ephemeral effects described above. The model can account for the fact that the ice maximum is not achieved on the same day of the ice cycle each year. It also recognizes that the length of the ice cycle will vary and that the timing of advance and retreat of the ice varies from year to year. This means that the annual cycle is not constrained to a fixed cyclical pattern, rather it is a pattern that allows both temporal dilation and contraction as well as amplitude modulation.

To show the utility of the model, we develop several different annual cycles, including one that is invariant, one that is adjusted for phase only and one that is adjusted for amplitude only. From the modeled annual cycles we define and extract the variability at the timescales mentioned. We conclude with a decomposition of the variability of SIE during 2016, the year of anomalous decay of SIE. The data are described in Sect. 2, and the model is defined and developed in Sect. 3. The results are presented and discussed in Sect. 3, while conclusions are made in Sect. 4.

2 Data

We used the Bootstrap Version 3 concentration fields (Comiso, 2017) from the NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration, Version 3 (Peng et al., 2013; Meier et al., 2017). These data were generated using the Advanced Microwave Scanning Radiometer – Earth Observing System (AMSR-E) Bootstrap Algorithm with daily varying tie points. They span the period 26 October 1978 to 31 December 2018 and are daily except prior to July 1987 when they are given every other day. Data are gridded on the SSM/I polar stereographic grid (25 km× 25 km). In addition to the alternate day observations from 1978 to 1987, there are a number of days and segments of days with no observations. In particular, there are no data between early December 1987 and mid-January 1988. Our methods do not require a complete temporal data record and naturally deal with missing data. As such we do not impute the missing days. The SIE used in our analysis was calculated using the conventional limit of the 15 % SIC isoline. Every grid poleward of the 15 % isoline is considered to be completely covered with ice.

Figure 1 shows the recorded total SIE (in grey dots) for each year from 1979 to 2017 and a smoothed representation of the traditional daily annual cycle (red). In this figure, day 0 on the horizontal axis represents the typical lowest SIE for the year, Julian day 50. We employ this convention for all of the time series figures used in this paper. The plot nicely illustrates the variation in the SIE from day-to-day and year to year.

Figure 1 Recorded sea ice extent (SIE) (grey) for each year compared to a smooth annual cycle (red) over a 365 d period. The horizontal axis is the day of the cycle and the vertical axis is sea ice extent in millions of square kilometers.

3 Methods and results: a statistical decomposition of sea ice extent

3.1 Annual cycle definition

In this section we give five ways to define an annual cycle in the sea ice extent. We start with the traditional definition of the annual cycle and progressively define annual cycles that are more sophisticated and that can represent more of the variation in the SIE over time. The second is an invariant annual cycle that retains the 365 d period of the traditional but incorporates the smooth functional form we might expect. The third adds amplitude variation to the invariant annual cycle so that the cycle itself varies from year to year with the amplitude of the year. The fourth adds phase variation to the invariant annual cycle, allowing it to capture the timing of the ice advance and retreat over each year. Finally, the fifth adds both amplitude and phase variation to the invariant annual cycle, allowing it to represent variation over time in both the amplitude and phase of the SIE.

Traditional annual cycle

Our decomposition of the sea ice extent starts with the traditional representation based on the annual cycle is as follows:

$$\begin{array}{}\text{(1)}& \text{extent}\left(t\right)=a\left[\text{doy}\right(t\left)\right]+\mathit{\alpha }\left(t\right)\text{where}t=T_{\mathrm{0}},\mathrm{\dots },T,\end{array}$$

where extent(t) is the extent on day t expressed as a decimal year (e.g., 1 February 2010 as 2010.08767) and doy(t) is the day of the year for t (e.g., 32). Most importantly, a is an annual cycle shape function with a(s) giving the annual cycle shape value for day of the year s. In this context, α(t) is the anomaly of the extent from the annual cycle on day t, T₀ is the first observed time and T is the last observed time. For the data in this paper, T₀=1978.833 and T=2019.000.

Within this representation, the annual cycle is traditionally estimated by a_T[s]:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{a}_T\left[s\right]={\displaystyle \frac{\mathrm{1}}{{\sum }_{t:\text{doy}\left(t\right)=s}\mathrm{1}}}\sum _{t:\text{doy}\left(t\right)=s}\text{extent}\left(t\right),\\ \text{(2)}& {\displaystyle }& {\displaystyle }\text{where}t=T_{\mathrm{0}},\mathrm{\dots },T,\end{array}$$

where ${\sum }_{t:\text{doy}\left(t\right)=s}\mathrm{1}=\mathrm{40}$ is the number of years of data.

This traditional estimate (a_T[s],) has a number of statistical issues that reduce its utility for examining the sea ice variability. Firstly, it is typically based on data for a subset of the satellite era (e.g., from 1979 onward). Currently, this is about 40 years of data, inducing intrinsic statistical variability into a_T[s] as an estimate of a[s]. This could be reduced by increasing the temporal range backward, by, for example, including data from the earlier satellite record (NIMBUS-5). Another option is to include information from proxy sources. However, this requires a large and sophisticated model-based reconstruction and we do not consider such methods further in this paper. Secondly, a_T[s] is computed separately for each day, ignoring the surrounding days. There is information in the temporally close days in the intuitive sense that days close to s, e.g., s−1 and s+1, will have similar, albeit not exactly the same, values. This information is ignored by a_T[s]. Thirdly, we expect a[s] to be smooth as a function of s so that changes in a_T[s] with s will be similar for days that are close. Fourthly, we expect that a_T[s] will “over fit” to the record, making the estimated anomalies from it smaller than the true anomaly, α(t), and that the annual cycle estimates will be more variable than the true annual cycle. This last issue is induced by the finite record and the estimates of the anomaly $\widehat{\mathit{\alpha }}\left(t\right)=\text{extent}\left(t\right)-a_T\left[\text{doy}\right(t\left)\right]$ will be statistically different than those of α(t). In summary, the traditional estimate, a_T[s], uses limited information, ignores other days, is not as smooth as we expect due to day-to-day variation, and over fits to the record.

Invariant annual cycle

It is possible that smoothing the data could be a solution to the statistical issues that arise from the way in which the traditional annual cycle is calculated. To address this we define an invariant annual cycle, a_I[s], which models a[s] as a cyclic cubic spline function (Wegman and Wright, 1983) of s. Specifically, a[s] is modeled as a piecewise cubic polynomial that has a continuous second derivative, is continuous, has continuous first and second derivatives at T, and best fits the recorded (satellite-observed) extents while being smooth. The specific criterion for the last feature is to choose a_I[s] to minimize the penalized-square error (PSE):

$$\begin{array}{ll}{\displaystyle }{\text{PSE}}_{\mathit{\lambda }}\left(a\right)=& {\displaystyle }\sum _{t=T_{\mathrm{0}}}^T\mathit{\left\{}\text{extent}\right(t)-a[\text{doy}\left(t\right)]{\mathit{\}}}^{\mathrm{2}}\\ \text{(3)}& {\displaystyle }& {\displaystyle }+\mathit{\lambda }\underset{\mathrm{0}}{\overset{\mathrm{365}}{\int }}{a}^{\prime \prime }[s{]}^{\mathrm{2}}\mathrm{d}s\mathit{\lambda }>\mathrm{0},\end{array}$$

where a^′′[s] is the second derivative of a[s] and λ is a smoothing parameter, chosen to balance the closeness of fit to the recorded values (the first term) with the smoothness of a[s] (the second term). Hence, choosing the function a[s] that minimizes PSE_λ(a) provides a balanced representation of the annual cycle. It prioritizes smoothness of a[s] over the closeness of fit of a[s] to the recorded extents. Note that the traditional estimator, a_T[s], is the minimizer with λ=0, i.e., with no penalty for lack of smoothness. The choice of λ is subjective. In this work we choose to maximize the ability to predict unrecorded extents. Specifically, we use generalized cross-validation (GCV) (Craven and Wahba, 1978) to choose and the R package mgcv of Simon Wood for analysis (Wood, 2004, 2017). The annual cycle obtained in this way is the optimal smoothest annual cycle chosen to minimize the mean-squared error (MSE) of SIE. Any trends are removed, and there is no adjustment for phase or amplitude. Figure 2a compares the traditional annual cycle (plotted from Julian day 50 in 2016 to day Julian day 49 in 2017) with the recorded SIE and the invariant annual cycle. The visual improvement is modest but, as shown in Table 1, the invariant annual cycle represents a 28.7 % improvement in the MSE compared to the traditional cycle. Note that both annual cycles overestimate the SIE in the retreat phase of the ice for 2016, which is known to be an anomalous year.

Amplitude-adjusted annual cycle

The invariant annual cycle has the same motivation as the traditional annual cycle while being a clear statistical and conceptual improvement over the traditional cycle. However, we argue that since it is also fixed by day of year, it may be too restrictive since it, like the traditional cycle, disguises the contributions of both amplitude and phase to the annual cycle. To address this we define a complementary annual cycle that is deformed each year in two ways. The first is amplitude in the sense that the yearly maximum and minimum extents may vary, but the shape of the daily extent may be invariant. We enable the annual cycle to vary from year-to-year as a parameterized function of the annual cycle shape function. Specifically, we define the amplitude-adjusted annual cycle, a_A[s,y], to satisfy the following expression:

$$\begin{array}{ll}{\displaystyle }\text{extent}\left(t\right)=a_A[& {\displaystyle }\text{doy}\left(t\right),\text{min}\cdot \text{extent}\left(\text{year}\right(t\left)\right),\\ \text{(4)}& {\displaystyle }& {\displaystyle }\text{max}\cdot \text{extent}\left(\text{year}\right(t\left)\right)]+\mathit{\alpha }(t),\end{array}$$

where

$$\begin{array}{}\text{(5)}& a_A[s,min,max]=u_A\left[s\right](max-min)+min,\end{array}$$

and year(t) is the year for t (e.g., 2010), max ⋅ extent(y) is the scale parameter giving the maximum extent for year y and min ⋅ extent(y)) is the scale parameter giving the minimum extent for year y. Here u_A[s] is an invariant annual cycle for the standardized extent. It is defined in an analogous way to the invariant annual cycle as a smooth function. Specifically, u_A[s] as a cyclic cubic spline function of s is chosen to minimize the penalized-square error:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{PSE}}_{{\mathit{\lambda }}_A}\left(u\right)=\\ {\displaystyle }& {\displaystyle }\sum _{t=T_{\mathrm{0}}}^T{\left\{{\displaystyle \frac{\text{extent}\left(t\right)-\text{min}\cdot \text{extent}\left(\text{year}\right(t\left)\right)}{\begin{array}{c}\text{max}\cdot \text{extent}\left(\text{year}\right(t\left)\right)\\ -\text{min}\cdot \text{extent}\left(\text{year}\right(t\left)\right)\end{array}}}-u\left[s\right]\right\}}^{\mathrm{2}}\\ \text{(6)}& {\displaystyle }& {\displaystyle }+{\mathit{\lambda }}_A\underset{\mathrm{0}}{\overset{\mathrm{365}}{\int }}{u}^{\prime \prime }[s{]}^{\mathrm{2}}\mathrm{d}s{\mathit{\lambda }}_A>\mathrm{0},\end{array}$$

where λ_A is a smoothing parameter with the same role as λ_I.

Figure 2 Comparison of annual cycle estimates: (a) a traditional and invariant cycle and (b) a traditional and amplitude- and phase-adjusted cycle. The horizontal axis is the day of the cycle and the vertical axis is sea ice extent in millions of square kilometers.

This annual cycle gives a different decomposition of the extent to the invariant annual cycle as it captures variation due to amplitude variation. Specifically, adjusting for amplitude results in a 55.2 % improvement in the MSE compared to the traditional cycle (see Table 1). Note that this allocates that component of the variation in extent due to amplitude variation to the annual cycle rather than the residual term, α(t) (see Eq. 4). The magnitude of the change clearly underscores the importance of amplitude variations in the definition of the annual cycle.

Phase-adjusted annual cycle

Another component of the annual cycle that is important is the phase. This is the timing of the maximum and minimum extents. It is important because it determines the length of the annual cycle and influences its shape. We enable the annual cycle to vary from year-to-year as a parameterized function of the phase of the annual cycle shape function, defining the phase-adjusted annual cycle, a_P[s], as follows:

$$\begin{array}{}\text{(7)}& \text{extent}\left(t\right)=a_P\left[\text{phase}\right(t\left)\right]+\mathit{\alpha }\left(t\right),\end{array}$$

where phase(t) is the phase-adjusted day of the year for t (e.g., 164). It is a smooth function of time that tells us what day of an invariant 365 d cycle the date t is. The function phase(t) is modeled here as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{phase}\left(t\right)=\mathrm{365}\\ \text{(8)}& {\displaystyle }& {\displaystyle }\times \text{Beta}\left({\displaystyle \frac{t-\text{min}\cdot \text{extent}\cdot \text{day}\left(\text{year}\right(t\left)\right)}{\begin{array}{c}\text{max}\cdot \text{extent}\cdot \text{day}\left(\text{year}\right(t)+\mathrm{1})\\ -\text{min}\cdot \text{extent}\cdot \text{day}\left(\text{year}\right(t\left)\right)\end{array}}};\mathit{\beta }\left(\text{year}\right(t\left)\right)\right)\text{(9)}& {\displaystyle }& {\displaystyle }\mathrm{\text{min}}\cdot \text{extent}\cdot \text{day}\left(\text{year}\right(t\left)\right)\le t\le \text{max}\cdot \text{extent}\cdot \text{day}\left(\text{year}\right(t\left)\right),\end{array}$$

where $\text{max}\cdot \text{extent}\cdot \text{day}\left(y\right)$ is the day of the year giving the maximum extent for year y and $\text{min}\cdot \text{extent}\cdot \text{day}\left(y\right))$ is the day of the year giving the maximum extent for year y. Here $\text{Beta}(p;\mathit{\beta }),\mathrm{0}\le p\le \mathrm{1},$ is the cumulative distribution function of a Beta(β) random variable parameterized by $\mathit{\beta }=({\mathit{\beta }}_{\mathrm{1}}>\mathrm{0},{\mathit{\beta }}_{\mathrm{2}}>\mathrm{0})$ and β(y) is the parameter value specific to year y.

Here a_P[s] is an invariant annual cycle for the extent (typically differing from a_I[s]). It is defined in an analogous way to the other invariant annual cycles as a cyclic cubic spline function of s chosen to minimize the penalized-square error:

$$\begin{array}{ll}{\displaystyle }{\text{PSE}}_{{\mathit{\lambda }}_P,\mathit{\beta }}\left(u\right)& {\displaystyle }=\sum _{t=T_{\mathrm{0}}}^T\mathit{\left\{}\text{extent}\right(t)-u[\text{phase}(t;\mathit{\beta }(\text{year}\left(t\right)\left)\right]{\mathit{\}}}^{\mathrm{2}}\\ \text{(10)}& {\displaystyle }& {\displaystyle }+{\mathit{\lambda }}_P\underset{\mathrm{0}}{\overset{\mathrm{365}}{\int }}{u}^{\prime \prime }[s{]}^{\mathrm{2}}\mathrm{d}s{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_P>\mathrm{0},\mathit{\beta }\left(\text{year}\right(t\left)\right)>\mathrm{0},\end{array}$$

where λ_P is a smoothing parameter, chosen to balance the closeness of fit to the recorded values (the first term) with the smoothness of u[s] (the second term). The minimization is also over the parameters $\mathit{\left\{}{\mathit{\beta }}_{\mathrm{1}}\right(y)>\mathrm{0},{\mathit{\beta }}_{\mathrm{2}}(y)>\mathrm{0}{\mathit{\}}}_{y=\mathrm{1978}}^{\mathrm{2018}}$.

The phase-adjusted annual cycle gives a different decomposition of the extent to the invariant annual cycle as it captures variation due to phase variation. It allocates the component of the variation in extent due to phase variation to the annual cycle rather than the residual term, α(t).

Surprisingly, the adjustment for phase shows even more improvement (63.9 %) in the MSE than that for the amplitude-adjusted annual cycle, indicating that the phase contributes more to the variability of the annual cycle of SIE than the amplitude. Most studies of Antarctic sea ice variability focus on the amplitude at maximum and minimum extents, but this analysis indicates that the phase (the timing of these extrema) is at least as important a contributor to the variability.

Amplitude- and phase-adjusted annual cycle

Finally, we can combine the amplitude and phase adjustment ideas to define an annual cycle that jointly adjusts for both. We define the amplitude- and phase-adjusted annual cycle (APAC), a_AP[s], as follows:

$$\begin{array}{ll}{\displaystyle }\text{extent}\left(t\right)=a_A[& {\displaystyle }\text{phase}\left(t\right),\text{min}\cdot \text{extent}\left(\text{year}\right(t\left)\right),\\ \text{(11)}& {\displaystyle }& {\displaystyle }\text{max}\cdot \text{extent}\left(\text{year}\right(t\left)\right)]+\mathit{\alpha }(t),\end{array}$$

where a_A and phase(t) are defined as in Eqs. (5) and (9). Note that they will be different functions as they are now jointly specified. As before, a_A[s] is modeled as a cyclic cubic spline function of s chosen to minimize the penalized-square error:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{PSE}}_{{\mathit{\lambda }}_{\text{APAC}},\mathit{\beta }}\left(u\right)=\\ {\displaystyle }& {\displaystyle }\sum _{t=T_{\mathrm{0}}}^T\mathit{\{}{\displaystyle \frac{\text{extent}\left(t\right)-\text{min}\cdot \text{extent}\left(\text{year}\right(t\left)\right)}{\text{max}\cdot \text{extent}\left(\text{year}\right(t\left)\right)-\text{min}\cdot \text{extent}\left(\text{year}\right(t\left)\right)}}\\ {\displaystyle }& {\displaystyle }-u\left[\text{phase}\right(t;\mathit{\beta }\left(\text{year}\right(t\left)\right)]{\mathit{\}}}^{\mathrm{2}}\\ \text{(12)}& {\displaystyle }& {\displaystyle }+{\mathit{\lambda }}_A\underset{\mathrm{0}}{\overset{\mathrm{365}}{\int }}{u}^{\prime \prime }[s{]}^{\mathrm{2}}\mathrm{d}s{\mathit{\lambda }}_{\text{APAC}}>\mathrm{0},\end{array}$$

where λ_APAC is a smoothing parameter. The minimization is also over the parameters $\mathit{\left\{}{\mathit{\beta }}_{\mathrm{1}}\right(y)>\mathrm{0},{\mathit{\beta }}_{\mathrm{2}}(y)>\mathrm{0}{\mathit{\}}}_{y=\mathrm{1978}}^{\mathrm{2018}}$. As for the other annual cycles (invariant, amplitude-adjusted, phase-adjusted), λ_APAC is chosen by generalized cross-validation.

Figure 2b compares the traditional annual cycle with the recorded SIE for 2016 and the APAC produced by this model for the same time period. The APAC is a much better fit to the recorded data and represents a large and significant improvement of 77.3 % in MSE (Table 1). Table 1 clearly demonstrates the value of having multiple successive definitions of the annual cycle when decomposing the variation in the daily annual cycle of SIE.

Table 1 Comparison of the various proposed annual cycles in terms of how well they explain the variation in daily SIE. Values are given as percentages of mean-squared error and the root-mean-squared error (RMSE).

The discussion above describes several different ways of defining the annual cycle of SIE. While an annual cycle adjusted for phase or amplitude only would not be the best estimate for the data, differences between them and the optimal estimated annual cycle (i.e., APAC) could reveal sources of variability in the daily SIE.

3.2 Analyzing variation: volatility, daily rate of change, anomalies and trend

Estimating the annual cycle using our model allows us to calculate statistics that reveal the underlying variability in the daily SIE. Below we decompose the sea ice variation on the time dimension, moving up the temporal scale from the very short term (the instantaneous variation) to the day-to-day variation, followed by the interannual variation and finally the multi-decadal variation, i.e., the trend.

The recorded sea ice extent will deviate from the true sea ice extent. This may be due to some combination of weather, artifacts of the satellite algorithm used for retrieval, and the electromagnetic spectrum across which the device or satellite is measuring, among other things. To represent this, we write the recorded SIE, SIE(t), as follows:

$$\begin{array}{ll}{\displaystyle }\text{SIE}\left(t\right)=& {\displaystyle }\text{extent}\left(t\right)+\mathit{ϵ}\left(t\right)\\ {\displaystyle }=& {\displaystyle }{a}_A\left[\text{phase}\right(t),\text{min}\cdot \text{extent}(\text{year}\left(t\right)),\\ \text{(13)}& {\displaystyle }& {\displaystyle }\mathrm{\text{max}}\cdot \text{extent}\left(\text{year}\right(t\left)\right)]+\mathit{\alpha }(t)+\mathit{ϵ}(t).\end{array}$$

The recorded SIE on any given day is then the sum of a number of components of variation – the annual cycle for that day, the yearly variation (anomaly) from the annual cycle and a residual term (*ϵ(t)). These are now discussed.

3.2.1 Volatility of the recorded sea ice extent

Here we introduce the term volatility to describe the instantaneous variation (or precision) in the recorded SIE as an approximation for the extent. Such variation may be due to ephemeral effects like those mentioned above.

Figure 3 Volatility of the recorded SIE for the NIMBUS-7 era (26 October 1978 to 20 August 1987) and the DMSP era (21 August 1987 to 2018). It is averaged over each day of the cycle in these eras. The units are given in millions of square kilometers. The purple curve is the day-to-day change in SIE from Fig. 4

Normally the standard deviation of the residual, ϵ(t) in Eq. (13), is represented as a constant over time. Here, however, we allow it to vary, explicitly representing it as a time-varying term or component. The volatility is therefore defined as the time series formed by the standard deviation of $\mathit{ϵ}\left(t\right),t=T_{\mathrm{0}},\mathrm{\dots },T$. It is a quantification of ephemeral effects. Effectively it shows the size and timing of the variability associated with factors like instrument error or noise in the recorded SIE.

To model the volatility, we specify a generalized autoregressive conditional heteroskedasticity (GARCH) model (Bollerslev, 1986) for the residual ϵ(t). The residual is split into a time-dependent standard deviation σ(t) representing the volatility and a series $z\left(t\right)\sim N(\mathrm{0},\mathrm{1})$:

$$\mathit{ϵ}\left(t\right)=\mathit{\sigma }\left(t\right)z\left(t\right).$$

Explicitly, the (squared) volatility is modeled as a weighted average of the past anomalies and (squared) volatilities:

$${\mathit{\sigma }}^{\mathrm{2}}\left(t\right)=\mathit{\omega }+\sum _{i=\mathrm{1}}^p{\mathit{\eta }}_i\mathit{ϵ}(t-i)+\sum _{i=\mathrm{1}}^q{\mathit{\psi }}_i{\mathit{\sigma }}^{\mathrm{2}}(t-i),$$

where the parameters η_i and ϕ_i represent dependency on the past residuals and volatilities, while the parameter ω represents a trend in volatility. The purpose of the dependency on past volatilities is to better represent periods of high or low volatility. We also specify an autoregressive moving average (ARMA) model for α(t) (Box and Jenkins, 1976; Hipel and McLeod, 1994) with ϵ(t) as the (time-dependent) error term. The model parameters were fit using maximum likelihood. The Bayesian information criterion (BIC) was used to select the model order (Ghalanos, 2019). The model orders were p=2 and q=2 (i.e., GARCH(2,2)) and auto-regressive moving average, ARMA(1,1), for the anomaly model. All models were fit using the R package rugarch (Ghalanos, 2019).

Figure 3 plots the average volatility in SIE, separating it into the two periods of time when different sensors were retrieving the data. There is some indication that the volatility is larger in the data recorded by the NIMBUS-7 Scanning Multi-channel Microwave Radiometer (SMMR) sensor (orange) than by the Defense Meteorological Satellite Program (DMSP) (black), especially at times of maximum SIE. This could be an effect of the sensor resolution (sensor footprint), which is actually smaller (higher resolution) in NIMBUS-7. These estimates adjust for the every-other-day sampling of the NIMBUS-7 sensor. Were this not adjusted for, the NIMBUS-7 values would be substantially higher than the DMSP. That said, there are some important similarities. Volatility is least at SIE minimum, larger at SIE maximum and largest late in the cycle when the ice is experiencing its largest rate of retreat. This latter characteristic is discussed below. The values from the DMSP era show that the volatility ranges approximately from 40 to 50 K km². These are relatively small values compared to the total SIE but are quite large compared to the typical grid cell size. The fact that the volatility is not constant over the cycle may be exploited to get a better understanding of contributors to overall variability in SIE.

3.2.2 Daily rate of change

It is useful to know the daily rate of change of SIE because it gives insight into the daily timing of growth (advance) and melt (retreat) of the sea ice. It is also an expression of the phase of the annual cycle. Contemporary trends in Antarctic sea ice are shown to be linked to the changes in the timing (phase) of advance and retreat (e.g., Stammerjohn et al., 2008). Note that the annual cycles have been defined as continuous in day. Hence, we can quantify the rate of change of total Antarctic SIE by the derivative of an annual cycle shape function, a[s]. The precise definition of the rate of change differs in the choice of annual cycle that is used. As an example, the rate of change for both the traditional and invariant annual cycles is plotted in Fig. 4, which shows the day-to-day changes in the SIE over the 365 d cycle. As might be expected, the overall pattern of the traditional (orange line) and invariant annual cycles (black line) are quite similar to each other. Both cycles show that the rates of growth and melt are variable over the cycle. However, compared to that of the invariant cycle, the day-to-day change in the traditional annual cycle is quite variable, making it difficult, if not impossible, to make precise statements about the timing of ice growth and decay. For example, around day 200 of the cycle, there is a reduction in the variability of the traditional annual cycle. This pause might be due to some idiosyncrasy in the data or it might be related to the relative stability of the ice extent in the region of the SIE maximum. The smooth monotonic day-to-day change of the invariant annual cycle shows that the day-to-day change is very close to zero, indicating that the latter reason is more likely. Therefore, the following comments are based on the day-to-day change in the invariant annual cycle.

The SIE minimum (day 0, Julian day 46) is coincident with the minimum growth rate. The ice advances, reaching the maximum growth rate by day 81 and maintaining this maximum growth rate for approximately 40 d before slowing to a minimum growth rate by day 225 (late September) of the cycle. Sea ice retreat begins at approximately day 225 and occurs quite rapidly compared to the advance, reaching a maximum rate at day 308 (late December) before slowing to a stop at day 365 (Julian day 46 or mid-February). The rates of advance and retreat of the ice are not constant over the annual cycle. The maximum rate of retreat of the ice is more than twice the maximum rate of advance. Figure 4 illustrates and more precisely defines a key characteristic of the Antarctic annual cycle, i.e., its asymmetry. The ice grows (advances) steadily over a much longer period than it decays (retreats). It has been suggested that this asymmetry in the annual cycle is a result of the influence of the semiannual oscillation (SAO) of the Antarctic circumpolar trough (Enomoto and Ohmura, 1990; Watkins and Simmonds, 1999) and an open-water (ice)–albedo feedback, with the latter being the main driver for the rapid retreat of sea ice (Ohshima and Nihashi, 2005). Ice budget analysis studies (Holland and Kwok, 2012; Holland, 2014; Holland and Kimura, 2016) indicate that surface winds are also important in the advance and retreat of the ice as they drive advection of ice and divergence within the pack. Recent modeling studies (Kusahara et al., 2019) suggest that ice advance is due chiefly to thermodynamic processes (except in the Ross Sea) while ice retreat is largely wind driven (or dynamic). Our study provides more precise information on the timing of advance and retreat and on the length of the two major stages of the ice cycle (ice growth: 225 d; ice retreat: 140 d) than can be obtained from monthly averaged data. This is significant because much of the variation in contemporary Antarctic SIE has been occurring at sub-monthly scales.

Figure 4 Day-to-day change in the annual cycle of sea ice extent for the traditional (orange) and invariant (black) annual cycles. The horizontal axis is the day of the cycle, and the vertical axis is change in sea ice extent in millions of square kilometers.

Taken together, the daily rate of change and the volatility (Figs. 3 and 4) show, (1) The timing of lowest volatility may be related to the fact that there is relatively little ice at minimum; (2) During the period when ice is advancing most swiftly, the volatility is low, responding to constant large-scale forcing; (3) During the period of slowing growth and maximum extent, volatility is high, perhaps due to the more frequent occurrence of storms during winter (Simmonds and Keay, 2000) causing fluctuations at the sea ice edge rather than within the pack where the sea ice concentration is at or close to 100 %. This effect of the storms may be magnified because at the ice maximum, the perimeter of the ice cover is also at or near its maximum, potentially allowing more ice area to be affected; (4) Volatility begins to decrease as the sea ice retreats; but (5) increases to its maximum value when the rate of retreat is largest. The late peak in volatility may be due to the dynamic nature of the retreat. Anecdotally, the sea ice extent anomalies of note tend to occur during the sea ice maximum and the period immediately following (Turner et al., 2017; Schlosser et al., 2018). The statistics examined here are suggesting that these anomalies are probably associated more strongly with dynamic forcing than thermodynamic.

3.2.3 Anomalies

The detection and analysis of anomalies (deviations from the annual cycle) is essential to the understanding of contributors to variability. Here we discuss three different but related types of anomalies. First there is the true anomaly, represented by α(t) in Eqs. (1), (11) and (13). This is the difference between the true SIE and the annual cycle, however it is defined. The true anomaly is the preferred anomaly but is unobtainable because of imprecision in measuring and retrieving the sea ice data. Second there is the raw anomaly, i.e., the difference between the observed (recorded) SIE and the annual cycle. Here we focus on a statistical estimate of the true anomaly, α(t), which we denote as $\widehat{\mathit{\alpha }}\left(t\right)$. The estimate is preferable to the raw anomaly as it adjusts for the volatility and should be closer to the true anomaly than the raw anomaly.

Figure 5 Comparison of anomalies from three annual cycle estimates for 2014–2018: the raw anomaly from the traditional annual cycle (black), the estimated anomaly from the invariant annual cycle (blue), and the estimated anomalies from the amplitude- and phase-adjusted annual cycle (red). The vertical axis is the anomaly in millions of square kilometers.

We estimate the true anomaly by using Eq. (13), rewriting it as follows:

$$\begin{array}{ll}{\displaystyle }\widehat{\mathit{\alpha }}\left(t\right)=\text{SIE}\left(t\right)-{\widehat{a}}_A[& {\displaystyle }\widehat{\text{phase}}\left(t\right),\text{min}\cdot \text{extent}\left(\text{year}\right(t\left)\right),\\ \text{(14)}& {\displaystyle }& {\displaystyle }\text{max}\cdot \text{extent}\left(\text{year}\right(t\left)\right)]-\widehat{\mathit{ϵ}}(t).\end{array}$$

We use the estimate of the APAC and compute $\widehat{\mathit{ϵ}}\left(t\right)$ from the GARCH model for the residual ϵ(t) from Sect. 3.2.1. The estimated anomaly is quite close to the recorded anomaly as $\widehat{\mathit{ϵ}}\left(t\right)$ is small in magnitude (see Figs. 3 and 7).

Figure 5 plots three types of anomalies: the raw anomaly from the traditional annual cycle and the estimated anomalies from the invariant and APAC. These show the last 5 years of the 42 years of satellite-observed data, 2014–2018. The anomalies of the three annual cycles are similar in sign; however, those for the APAC tend to be smaller. The similarity in sign is expected, and the smaller size of the APAC anomalies arises because the APAC is a much better fit to the recorded data. The anomalies for the traditional and invariant annual cycles are not significantly different from each other in size. This is expected given the small difference between the two shown in Fig. 2. We can clearly see the large negative anomaly in SIE at the end of 2016. The negative anomaly is larger in the traditional and invariant annual cycles than in the APAC, demonstrating that the APAC is a better fit to the recorded SIE; therefore, the anomaly is expected to be smaller.

3.2.4 Trend

The trends in SIE for both the Arctic and Antarctic have been the subject of much study. Most studies assume a linear trend and employ a linear model of the monthly data to estimate those trends (e.g., Parkinson and Cavalieri, 2012). Instead, we remove this assumption of linearity and model the trend in the daily data as a thin plate regression spline function of time (Wood, 2003). We added a term to our model for the SIE representing this curvilinear trend and jointly estimate it by minimizing the PSE (penalized-square error):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{PSE}}_{{\mathit{\lambda }}_I,{\mathit{\lambda }}_{\text{trend}}}(a,\text{trend})=\\ {\displaystyle }& {\displaystyle }\sum _{t=T_{\mathrm{0}}}^T\mathit{\left\{}\text{extent}\right(t)-\text{trend}(t)-a[\text{doy}\left(t\right)\left)\right]{\mathit{\}}}^{\mathrm{2}}+{\mathit{\lambda }}_I\underset{\mathrm{0}}{\overset{\mathrm{365}}{\int }}{a}^{\prime \prime }[s{]}^{\mathrm{2}}\mathrm{d}s\\ \text{(15)}& {\displaystyle }& {\displaystyle }+{\mathit{\lambda }}_{\text{trend}}\underset{T_{\mathrm{0}}}{\overset{T}{\int }}{\text{trend}}^{\prime \prime }(t{)}^{\mathrm{2}}\mathrm{d}t{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_I>\mathrm{0},{\mathit{\lambda }}_{\text{trend}}>\mathrm{0},\end{array}$$

where trend^′′(t) is the second derivative of trend(t) at time t and λ_trend is a smoothing parameter specific to the trend and is chosen to balance the closeness of fit to the recorded values using generalized cross-validation (Wood, 2004). The last term also captures the beginning and end time smoothing.

Figure 6 Three estimates of the trend in the recorded SIE represented in terms of the amount of SIE associated with the change. The blue line is the linear trend estimated for data from 1 January 1979 to 31 December 2017. The red line is the linear trend estimated for data from 1 January 1979 to 31 December 2018. The black line is the curvilinear trend estimated for data from 26 October 1978 to 31 December 2019. The dashed lines are the 95 % pointwise confidence bands for the smooth curvilinear trend.

The curvilinear trend in SIE for 1979–2015 and 1979–2018 derived using this method is illustrated in Fig. 6 along with the linear estimates of the trend. The latter assumes the same model as Eq. (15) except it constrains trend(t) to be linear. While there is a small positive linear trend, as has been reported in the literature (e.g., Parkinson and Cavalieri, 2012; Turner et al., 2015a), Fig. 6 shows that there is strong nonlinearity in the trend. There are strong decadal differences. For example, in the 1980s the trend was largely negative, while from 1990 to the mid-2000s there were a number of short-term fluctuations with opposing signs. It seems clear from Fig. 6 that the reported positive trend in total Antarctic SIE is due largely to the positive trend that began at the end of the first decade of the 21st century and continued until 2016. The anomalously low SIE experienced since 2016 had the effect of reducing the slope of the linear trend by almost 50 % from 13 860 to 6068 km² yr⁻¹. If the trends were linear they would be statistically significantly positive. The nonlinearity of the daily SIE trend in this analysis is consistent with that discussed by Simpkins et al. (2013) in their analysis of changes in the magnitudes of the sea ice trends in the Ross and Bellingshausen Seas. We note also that use of the daily data adjusted for amplitude and phase potentially allows a better estimate of the trend than monthly averaged values.

Figure 7 (a) Decomposition of the sea ice extent during 2016 into various components of its variation, including separate amplitude and phase components. (b) The invariant and amplitude-adjusted annual cycles for 2016. Day of cycle is on the horizontal axis: day 0 is Julian day 52. Anomalous SIE in millions of kilometers squared is on the vertical axis.

Even within the context of nonlinearity, the anomalously low SIE represents a dramatic negative adjustment to Antarctic SIE (Schlosser et al., 2018; Parkinson, 2019), prompting questions about whether or not this represents a change in state instead of a fluctuation due to natural variability. The current length of record does not allow much more than speculation. However, we can decompose the annual cycle of 2016 into the various components of variation that we have identified in this paper. This is illustrated in Fig. 7. The daily values of the components are plotted against the anomaly in SIE, showing how much they contributed to the SIE anomaly. The decomposition is sequential with the amplitude component extracted before the phase component.

The decomposition shows that the curvilinear trend (green) for 2016 is small and positive early in the cycle and becomes strongly negative later in the year, making a large contribution to the negative anomaly during this time of rapid change, as identified in Fig. 6. The raw anomaly (black), the difference between the recorded SIE and the APAC (the anomaly that includes the volatility) and its smoothed version, the estimated anomaly (orange), is small and did not make a consistent contribution to the anomaly in SIE over the year. Smoothing removed the “noise”, which might be due to instrumentation and leaving behind a truer variation between the recorded SIE and the expected SIE (i.e the APAC). The amplitude (blue) made a small but consistently negative contribution to the anomalously low SIE in 2016. Interestingly, the main contributor to the anomalous SIE was the phase. That is, the phase contributed to a small positive anomaly during the growth stage of the cycle (the growth was slightly ahead of phase) and a strongly negative anomaly during retreat, indicating that the timing of retreat of the ice was earlier than normal and the ice retreated faster than normal. The sum of these components (including the invariant annual cycle (pale blue) is the recorded SIE for 2016. The decomposition shows that the difference between the recorded SIE and the traditional and invariant cycles seen in Fig. 2a is mostly due to phase, as is also seen in Fig. 7b.

4 Conclusions

Variability in the annual cycle of Antarctic sea ice extent is dominated by the amplitude and phase of the cycle. In this study, we examined the variability in the annual cycle of total Antarctic sea ice extent (SIE) in detail at timescales ranging from instantaneous to day-to-day, interannual, and multi-decadal trends, thus offering a complete picture of the temporal variation in the sea ice. To facilitate this analysis, we developed first a statistical and mathematical model of the annual cycle in which the amplitude and phase, the two major contributors to its variability, are allowed to vary. This is contrary to traditional methods that restrict the variation in amplitude and phase, thus limiting their contribution to the variability. We define a number of complementary annual cycles – the invariant, which is an optimally smoothed annual cycle with no adjustments for phase or amplitude, an annual cycle that adjusted for phase only, another adjusted for amplitude only, and one that is adjusted for phase and amplitude (APAC). Each of these annual cycles represent clear conceptual and statistical improvements over the traditional method of calculating the annual cycle, with the APAC showing the most improvement. We propose the APAC as a substitute for the traditional method. However, the differences between the other annual cycles and the APAC reveal sources of variability in the daily SIE. For example, comparing the annual cycles adjusted for phase only and amplitude only revealed that the phase contributes more to the variability in the annual cycle than the amplitude.

The timescales into which the variability of SIE was decomposed allow useful interpretations of the factors that give rise to the variability. Using the volatility defined and described here for the first time, we show how much of the total SIE is due to transient effects (such as satellite effects and algorithmic artifacts). We also show how those effects vary over the annual cycle, and in the process we note that there are differences in the volatility (and hence uncertainty) that arise because of sensor type. The daily rate of change in SIE allows a precise definition of the timing and rate of advance and retreat of the sea ice, a quality that is very important given that much of the contemporary variability in Antarctic sea ice occurs at sub-monthly scales. Combination of the information given by the volatility and daily rate of change suggests that the volatility is lowest when the sea ice is at its minimum and highest during the time of the maximum rate of retreat. Given that the rapid rate of retreat of the ice has been associated with dynamic processes, this suggests that the peak in volatility at the end of the cycle is due to ephemeral effects associated with dynamic forcing.

To look at the interannual timescale, we defined and estimated several different but related anomalies, i.e., measures of deviation from the annual cycle, that may be used to evaluate the contributions to Antarctic sea ice variability from sources (local, oceanic and atmospheric) other than the large-scale sources that control cyclical, amplitude and phase changes. These show that our proposed annual cycle, the APAC, is a better fit to the recorded SIE.

We established that the trend in daily total Antarctic SIE over time is strongly nonlinear and that the linear estimates are weak and dependent on a positive trend that began in 2011 and ended in 2016. Interestingly, our decomposition of the annual cycle of 2016 into the components of variation defined in this paper shows that the main contributor to the anomalous SIE was the phase. That is, the anomalously low SIE was due mainly to the fact that the retreat began earlier than normal and was faster than normal. The amplitude made a much smaller negative contribution that did not vary much over the year.

We used the daily total Antarctic SIE in this analysis. However, sea ice variability around Antarctica is strongly regional, and the annual cycles of these regions are markedly different from each other and are also changing. The model-estimated annual cycles and the timescale decomposition presented here will facilitate examination of the regional variability of Antarctic sea ice. Finally, although our method was developed on Antarctic SIE, this decomposition methodology is applicable to a wide range of climatic variables (e.g., temperature, Arctic sea ice extent) that experience an annual cycle.