A new tracking algorithm for sea ice age distribution estimation

A new algorithm for estimating sea ice age (SIA) distribution based on the Eulerian advection scheme is presented. The advection scheme accounts for the observed divergence or convergence and freezing or melting of sea ice and predicts consequent generation or loss of new ice. The algorithm uses daily gridded sea ice drift and sea ice concentration products from the Ocean and Sea Ice Satellite Application Facility. The major advantage of the new algorithm is the ability to generate individual ice age fractions in each pixel of the output product or, in other words, to provide a frequency distribution of the ice age allowing to apply mean, median, weighted average or other statistical measures. Comparison with the National Snow and Ice Data Center SIA product revealed several improvements of the new SIA maps and time series. First, the application of the Eulerian scheme provides smooth distribution of the ice age parameters and prevents product undersampling which may occur when a Lagrangian tracking approach is used. Second, utilization of the new sea ice drift product void of artifacts from EUMETSAT OSI SAF resulted in more accurate and reliable spatial distribution of ice age fractions. Third, constraining SIA computations by the observed sea ice concentration expectedly led to considerable reduction of multi-year ice (MYI) fractions. MYI concentration is computed as a sum of all MYI fractions and compares well to the MYI products based on passive and active microwave and SAR products.

1 Introduction

Sea ice age (SIA) is one of the components of the essential climate variable (ECV) for sea ice as defined by the Global Climate Observing System (GCOS) (WMO, 2016). It is an important climate change indicator in polar regions, which describes the sea ice cover state in addition to its concentration and thickness. More generally, the changes in SIA can be seen as a proxy of sea ice decline in Arctic.

Numerous studies have been focused on the estimation and evolution of the Arctic SIA over the last decades (Fowler et al., 2004; Maslanik et al., 2007, 2011; Rigor and Wallace, 2004; Tschudi et al., 2010, 2016c). They all report that the amount of old/thick ice in the Arctic Ocean has been decreasing dramatically over the last decades. For example, the very old ice (> 4 years old) comprised 20 % of the ice pack in 1985, decreasing to 3 % in 2015, while the 3- to 4-year-old ice decreased by one-third over the same time period (Perovich et al., 2015). Kwok (2004), Polyakov et al. (2012) and Kwok and Cunningham (2015) reported significant decline in the area covered by multi-year ice (MYI) in the central Arctic since 1999, to which they associated a significant decrease in the mean sea ice thickness and volume in the same region.

Sea ice age is an important parameter of the Arctic ocean system and may be a good indicator of significant changes in the dynamical and thermodynamical regimes that have taken place since the beginning of the current century, such as sea ice thinning and faster ice drift. The younger seasonal ice cover present in the Arctic nowadays is in general thinner, which makes it more vulnerable to break up, deformation and drift under the actions of waves, winds and currents. This has been confirmed by calculating the trends in sea ice drift and deformation over the last 3 decades, which are significant and respectively about 13 and 50 % per decade (Rampal et al., 2009). If more mobile, sea ice might also be exported more easily and rapidly out of the Arctic basin, e.g., through Fram Strait, especially the MY sea ice pack located north of Greenland and Ellesmere Island. A more fragmented sea ice cover is also more vulnerable to summer melt. For example, during the winter of 2007–2008 large chunks of MYI resulting from the intense break up of the thickest part of the Arctic ice cover located in the north of the Canadian Archipelago were melted away while drifting in the Beaufort Sea within the few following summer months (Kwok and Cunningham, 2010). A positive feedback on the loss of MYI due to the more frequent sea ice cover fracturing, increased sea ice drift and enhanced melting rate is therefore potentially activated, contributing to the observed negative trend in SIA.

It is the purpose of the present paper to describe a method and a derived dataset that allow us to shed more light on the development of the age distribution of the Arctic sea ice. For this purpose, we have taken advantage of some new datasets on sea ice drift and concentration developed and distributed by the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI SAF). In addition, we have developed a new Eulerian scheme of advection supported by the Sea Ice Climate Change Initiative (SICCI) project of the European Space Agency (ESA). These improvements have allowed us to avoid the problem of the tracers dispersion (Rybak and Huybrechts, 2003) and to produce a new SIA dataset which in each grid box contains not only the age of the oldest ice, but the actual age distribution provided as fractions of ice of different age categories (hereafter referred to as “sea ice age fractions”). The dataset will be presented and compared with earlier attempts to map Arctic SIA as well as with the standard products for sea ice type classification from scatterometer and microwave radiometer observations.

2 Data

2.1 Sea ice drift

Information on sea ice motion was acquired from two sources. First, the National Snow and Ice Data Center (NSIDC) sea ice drift (SID) product v.0116 (Tschudi et al., 2016c) was downloaded from the NSIDC portal (Tschudi et al., 2016b). Weekly SID fields from NSIDC were accessed on an Equal Area Scalable Earth (EASE) grid (Brodzik et al., 2012) with 25 km spacing for the period from October 1978 to December 2015. As pointed out by Szanyi et al. (2016) this ice drift product contains artifacts due to the composition of the ice drift derived from satellite data, observed by in situ buoys and predicted using a simple free drift model. On short timescales, these artifacts result in large openings in the predicted distribution of MYI which are filled with the first-year ice (FYI). On longer timescales the openings get bigger and the shape of the predicted ice age distribution may be significantly distorted.

Figure 1 Sea ice age fractions for 1 January 1985 in the Greenland and Lincoln seas. It is clearly seen that in many areas the fraction of 8YI does not exceed 20 %, but the corresponding SIA map shows that the entire pixel is assigned to contain 8YI.

The second SID product was produced by the Ocean and Sea Ice Satellite Application Facility (OSI SAF) High Latitude Processing Center. The low-resolution sea ice drift product from the EUMETSAT OSI SAF is operational since 2009. It implements the continuous maximum cross-correlation algorithm of Lavergne et al. (2010) for retrieving ice drift vectors from daily composited maps of various medium resolution passive and active microwave (PAMW) satellite sensors such as the SSMIS, AMSR2 and ASCAT. A multi-sensor analysis is also provided (Lavergne, 2016a). The retrieval of SIA requires sea ice drift information during summer. This was recently achieved by the OSI SAF product (Lavergne, 2016b) using the 18.7 GHz channels of the GCOM-W1 AMSR2 instrument (Kwok, 2008). For this study, a short re-processed record of OSI SAF ice drift product (October 2012–May 2017) was accessed, complemented by the operational product available from http://osisaf.met.no (May 2017 onwards). The SID product from OSI SAF was filtered component-wise with 3 × 3 pixels median filter and then upscaled using linear interpolation onto a grid with 10 km spatial resolution. Gaps in the product were filled with nearest neighbor values.

2.2 Sea ice concentration

Sea ice concentration (SIC) product v.1.6 (Tonboe et al., 2017a) was produced by the OSI SAF portal at 1-day temporal resolution and 10 km spatial resolution for the period September 2012–September 2017 covering the Northern Hemisphere. Validation of the product indicates that SIC can be retrieved with 10 % accuracy with slight seasonal variations (Tonboe et al., 2017b).

2.3 Sea ice type (MYI concentration)

Sea ice types can be discriminated with PAMW satellite observations since the physical signatures of sea ice change significantly after the influence of summer melt and brine rejection. Therefore sea ice that has survived at least one summer melt is referred to as multi-year ice, and seasonal ice is referred to as first-year ice.

The algorithm of Environment Canada Ice Concentration Extractor (ECICE) (Shokr et al., 2008) can combine several different observation inputs, e.g., passive and active microwave satellite data. Combined PAMW data can help to identify MYI; however, the retrieval shows flaws under specific weather conditions. The improved MYI concentration retrievals developed at University of Bremen is based on the ECICE algorithm using brightness temperatures from AMSR2 and radar backscatter from C-band scatterometer ASCAT. Two correction schemes were applied to the initial MYI concentration: one using temperature records from atmospheric reanalysis to identify MYI anomalies caused by warm spells and replace with interpolated MYI concentrations (Ye et al., 2016a) and the other utilizing mainly ice drift records to constrain the MYI changes within a plausible contour (Ye et al., 2016b).

The improved MYI concentrations were provided as gridded products on polar stereographic grid with 12.5 km spatial resolution for the winter months (November–April) of 2013–2017. Compared to the MYI dataset, for which the corrections were developed, the retrievals from AMSR2 and ASCAT are much coarser (12.5 km vs. 4.45 km). The coarser resolution requires adjustment of the thresholds in the drift correction, which results in suboptimal performance of the drift correction and therefore could lead to unexpected problems in the final dataset.

The OSI SAF sea ice type is another algorithm that combines both passive and active microwave data in a Bayesian approach (Aaboe et al., 2017) that computes the probability of occurrence of the most likely ice type – first-year ice or multi-year ice. The OSI SAF type product is a near-real-time product that has been operational since 2005 with data available from the OSI SAF High Latitude processing center (http://osisaf.met.no). Since the start of operational production the algorithm has been upgraded several times by including new sensors and improving the methods. For the period of interest in this study, 2013–2017 summer, the retrieval is based on scatterometer data from the ASCAT instrument and atmospherically corrected brightness temperatures from SSMIS. From October 2015 an algorithm upgrade resulted in the implementation of dynamical training data updated on a daily basis in order to replace the previous fixed training data. The ice type product is provided for the winter period of October until mid-May on polar stereographic projection with 10 km grid resolution.

2.4 Sea ice age

Information on SIA independent of the method introduced here was acquired from the NSIDC portal (Tschudi et al., 2016a). Weekly SIA fields were provided at EASE grid (Brodzik et al., 2012) with 12.5 km spatial resolution for the period from October 1978 to December 2015. The gridded SIA product is generated from positions of the virtual Lagrangian ice parcels (initially released on a regular grid) (Fowler et al., 2004) when the grid cell is assigned the age of the oldest parcel. This does not take into account fractions of younger ice present in the same cell. It may happen that most of the drifters within a cell are very young and only one drifter is old but nonetheless the entire cell is still assigned to be the oldest ice. Such a case is shown in Fig. 1 for 1 January 1985. For generating this figure we have implemented the NSIDC ice age algorithm and applied it to the NSIDC ice drift product starting from 1978. Concentration of each ice age category is calculated as a relative number of ice parcels with respective age.

3 Methodology

The SIA computation is implemented in two stages as explained in detail below. First, SIA fractions are independently advected using the satellite-derived sea ice drift and concentration observations and the Eulerian advection scheme. Second, the SIA is computed, accounting for concentration of each ice age fraction.

3.1 Eulerian advection scheme

At each time step the observations of sea ice drift velocity components (U and V) and sea ice concentration (C_OBS) are provided as gridded products on the polar stereographic grid with the same spatial resolution with pixel size in x and y directions equal to R_x and R_y, correspondingly. Uncertainties in U, V and C_OBS products impact the accuracy of the end product but are not considered in the present study. We assume that a sample ice parcel in a pixel has shape of one grid cell and an initial ice concentration C and corresponding velocity components U and V. An example is presented in Fig. 2, where the ice parcel originating from the pixel 31 is shown as a blue square. The parcel drifts from point A over time Δt and the coordinates of the destination point B can be computed as follows:

$$\begin{array}{r}{\displaystyle }{X}_B=X_A+U\mathrm{\Delta }t\\ {\displaystyle }{Y}_B=Y_A+V\mathrm{\Delta }t.\end{array}$$

Figure 2 Scheme of Eulerian advection of sea ice. Points A and B denote start and end of drift of an ice parcel shown by the blue square.

The areal fluxes of sea ice out of the donor pixel 31 into recipient pixels 12, 13, 22 and 23 are calculated as follows:

$$\begin{array}{r}{\displaystyle }{F}_{\mathrm{31}\mathit{▸}\mathrm{12}}=C_{\mathrm{31}}ad,\\ {\displaystyle }{F}_{\mathrm{31}\mathit{▸}\mathrm{13}}=C_{\mathrm{31}}bd,\\ {\displaystyle }{F}_{\mathrm{31}\mathit{▸}\mathrm{22}}=C_{\mathrm{31}}ac,\\ {\displaystyle }{F}_{\mathrm{31}\mathit{▸}\mathrm{23}}=C_{\mathrm{31}}bc,\end{array}$$

where a, b, c and d are sides of the rectangles formed by intersection of the displaced donor pixel and recipient pixels (Fig. 2). The recipient pixels are selected based on position of point B – the four closest pixels become recipients. Concentration of sea ice at the next step in a recipient pixel (e.g., pixel 13 in Fig.  2) ($C_{\mathrm{13}}^{*}$) is calculated as

$$\begin{array}{}\text{(3)}& C_{\mathrm{13}}^{*}=F_{\mathrm{21}\mathit{▸}\mathrm{13}}+F_{\mathrm{22}\mathit{▸}\mathrm{13}}+F_{\mathrm{31}\mathit{▸}\mathrm{13}}+F_{\mathrm{32}\mathit{▸}\mathrm{13}}\end{array}$$

or in generalized form:

$$\begin{array}{ll}\text{(4)}& {\displaystyle }{C}_R^{*}=\sum _{\mathrm{4}\text{donors}}{C}_{\text{D}}& {\displaystyle }\left(\mathrm{1}-{\displaystyle \frac{|X_R-X_{\text{D}}-U\mathrm{\Delta }t|}{R_x}}\right){\displaystyle }& {\displaystyle }\left(\mathrm{1}-{\displaystyle \frac{|Y_R-Y_{\text{D}}-V\mathrm{\Delta }t|}{R_y}}\right),\end{array}$$

where C_D is the concentration of a donor pixel, X_D and Y_D are the coordinates of a donor pixel, X_R and Y_R are the coordinates of the recipient pixels and R_x and R_y are pixel size.

Figure 3 Example of advection of an ice fraction from 1 October 2012 to 1 October 2013 shown for every third month.

When the observed sea ice drift field diverges then a gap appears between the ice parcels (e.g., in pixel 23). If the sum of fluxes (Eq. 4) into a recipient pixel falls behind the actually observed concentration then this is interpreted as opening of leads, freezing of water in the leads and generation of new ice. Then, ice is presented by two fractions: the partial concentration of the older ice C_OI is computed as the sum of incoming fluxes (Eq. 4) and the partial concentration of the young ice C_YI is computed as a remainder:

$$\begin{array}{}\text{(5)}& C_{\text{YI}}=C_{\text{OBS}}-C_{\text{OI}}.\end{array}$$

When the observed sea ice drift field converges then the advected ice parcels overlap (e.g., in the pixel 12) and the sum of fluxes into a recipient pixel may exceed the concentration observed by satellites. This is interpreted as generation of pressure ridges and increase in thickness of the older fraction of sea ice. In that case the total ice concentration at the next step is assigned to be the observed ice concentration.

During the freeze-up period the observed concentration increases and becomes higher than the sum of fluxes into a recipient pixel. This is also interpreted as generation of young ice (YI) only and Eq. (5) is used for computing YI concentration. During melting the observed concentration decreases and it may become less than the total predicted concentration (sum of OI and YI fractions). In that case the YI concentration is decreased first, and if YI is absent then the OI concentration is decreased.

3.2 Advection of sea ice age fractions

Advection of a SIA fraction is initiated on 10 September, the approximate date when ice extent reaches minimum in the Arctic, area of FYI is zero and all observed sea ice is MYI by definition. Initialization before this date has little impact, but initialization after this date increases risk of considering the observed FYI as MYI. The age of each ice fraction is increased by 1 year on 10 September of each consecutive year of advection.

In our study we initiated SIA calculation on 1 October 2012 when continuous high-quality observations of ice drift started to be available from AMSR2 (Fig. 3). For this date the total observed concentration was computed as the minimum concentration during the period 1 September–1 October 2002. Propagation of ice age fractions from later years started from 10 September or respective years.

We did not know the spatial distribution of ice of different age within the pack at the first moment of time (1 October 2012), but we can postulate that all observed ice is at least in the second-year ice (2YI) category. We have to make a simplification: the concentration of 2YI on 1 October 2012 is assigned to be equal to the total observed concentration: C_2YI=C_TOT (Fig. 3a). Understandingly, this simplification allows us to initialize the SIA algorithm but does not allow to provide full SIA frequency distribution before a sufficiently long spin-up period (e.g., 5 years). Such a spin-up period is also needed for the NSIDC product. Then the developed advection scheme is applied on a daily basis utilizing the available daily ice drift and concentration products (Fig. 3b–e). From 1 October 2012 to 10 September 2013 (Fig. 3a–d) the advected fraction contains 2YI but on 10 September 2013 the age of this fraction is increased by 1 year and later it becomes the fraction of the third-year ice (3YI) (Fig. 3e).

On 10 September 2013 the total ice concentration is the sum of all MYI fractions and therefore the initial concentration of the 2YI is calculated as

$$\begin{array}{}\text{(6)}& C_{\text{2YI}}=C_{\text{TOT}}-C_{\text{3YI}},\end{array}$$

where C_TOT is the ice concentration observed by satellites and C_3YI is the ice fraction advected from 1 October 2012. Both 2YI and 3YI fractions are advected further using the developed scheme and after N years of advection on each of 10 September the concentration of 2YI is calculated as

$$\begin{array}{}\text{(7)}& C_{\text{2YI}}=C_{\text{TOT}}-\sum _i{C}_{\text{i}},\end{array}$$

where C_i is the concentration of the ice fraction advected from previous years.

3.3 Computation of sea ice age

At any moment of time, sea ice in each cell is characterized by a vector of concentrations of ice fractions of various age with the sum equal to the total sea ice concentration (Fig. 4). In other words the new algorithm provides ice age probability distribution. A single number that characterizes a probability density function can be computed in several standard ways: minimum, maximum, mean, median, mode, etc. We propose that SIA can also be computed from the concentration of ice age fractions using several approaches depending on the preferred definition of the ice age (Fig. 5):

• Maximum age is the age of the oldest ice fraction that exceeds a given threshold (e.g., 5 % concentration).

• Average age is the mean of age of the ice fraction that exceeds a given threshold.

• Modal age is the age of the ice fraction that has the highest concentration.

• Median age is the age of the ice fraction that splits the age distribution into two equal parts. Linear interpolation between the ice age fractions or approximation of ice distribution by a predefined function may be needed to compute median age correctly. More generally, median (the value of the 50th percentile) can be replaced by any percentile.

• Weighted average age is the average of ages of individual ice fractions weighted by their concentrations:

  $$\begin{array}{}\text{(8)}& A={\displaystyle \frac{{\sum }_{\text{f}}{A}_{\text{f}}{C}_{\text{f}}}{{\sum }_{\text{f}}{C}_{\text{f}}}},\end{array}$$

  where A_f is age of a sea ice fraction and C_f is concentration of a sea ice fraction.

Figure 4 Maps of sea ice age fractions on 1 March 2015. (a) Observed total sea ice concentration; (b–d) fractions of advected multi-year ice; (e)  FYI fraction. Color bar represents the sea ice age fraction concentration.

Figure 5 Maps of sea ice age on 1 March 2015 computed with different methods. (a) Maximum ice age; (b) average ice age; (c) modal ice age; (d) weighted average ice age; (e) median (percentile 50) ice age; (f) NSIDC SIA. Color bar represents the sea ice age in years.

4 Results

The developed advection scheme (see Sect. 3.1) was applied to the sea ice drift and concentration products available at OSI SAF from 1 October 2012 to 1 March 2017. We compared our products to the already available sea ice products including SIA maps from NSIDC: time series and maps of MYI from NSIDC, Bremen University, OSI SAF and Denmark Technological University (DTU). The results of comparison presented below can be used for indirect assessment of our SIA algorithm.

4.1 Comparison with sea ice age product from NSIDC

In order to provide stepwise illustration of the SIA product improvements we have used four different combinations of forcings and algorithms to produce a SIA map for the 1 January 2016. First, we have implemented the NSIDC advection scheme (Lagrangian) and the SIA algorithm (age of the oldest parcel) and applied it to the ice drift product from NSIDC. Second, we have applied the NSIDC algorithms to the ice drift from OSI SAF. Third, we have applied our Eulerian advection scheme and SIA algorithm to the ice drift from OSI SAF but without accurately accounting for SIC. In this experiment all pixels with SIC below 15 % were assigned to be open water and other pixels contain 100 % ice. Finally, we have used the new advection scheme, the OSI SAF ice drift, and fully accounted for SIC.

Comparison of the generated ice age products is presented in Fig. 6. Panels a and b contain maps of maximum age computed with the NSIDC algorithm and panels c and d contain ice age estimated with the weighted average. The comparison indicates that spatial distributions of SIA are similar across four combinations of forcings and algorithms only in general. The change of the forcing significantly affects distribution of SIA (Fig. 6a and b) as was also illustrated in Szanyi et al. (2016). The belt of very old sea ice in the Lincoln Sea almost disappears, large areas with old ice in the central Arctic become more pronounced, long sleeves of 2YI with interspersing of older ice stretch across the central Arctic towards the Laptev Sea and along the Beaufort Sea coast into the Chukchi Sea, and vast areas of 2YI in the European Arctic become more homogeneous.

Figure 6 Comparison of SIA for the 1 January 2016 calculated with the following combinations of forcings and algorithms: (a) NSIDC drift and NSIDC algorithm; (b) OSI SAF drift and NSIDC algorithm; (c) OSI SAF drift and SICCI algorithm without SIC; (d) OSI SAF drift and SICCI algorithm with SIC.

Figure 7 Comparison of MYI concentration estimated from the NSIDC SIA product (a), from the SICCI SIA product (b), downloaded from OSI SAF (c) and from University of Bremen (d). Daily maps are generated for 31 December of the years from 2012 to 2016 and for 29 March 2017 from available data. Red contours denote extent of MYI in the SICCI product at 5 % concentration.

Figure 8 Seasonal dynamics of area of sea ice age fractions (filled plots) and multi-year ice area (dots) derived from the NSIDC product (a), from the SICCI product (b), from the UB product (black dots) and from the OSI SAF (yellow dots).

The change of the advection scheme does not affect the overall spatial distribution of MYI significantly (Fig. 6b and c) but generates a much smoother picture without discontinuities. Using the weighted average for SIA computation dramatically decreases the observed age (Fig. 6c). Few small spots where SIA reaches 4.5–5 years are observed only near the Canadian Archipelago and in the central Arctic. A stripe of high SIA appears along sea ice edge in the Fram Strait. This ice was advected from the central Arctic but fractions of younger ice have melted and only fractions of older ice remain, making the weighted average value high.

When the new algorithm takes SIC into account then SIA is decreasing even more and the ice older than 4 years is observed only near the Canadian Archipelago coast (Fig. 6d). Moreover, extent of MYI is also decreased and the sleeves of 2YI towards the Laptev and Chukchi seas almost disappear.

4.2 Intercomparison of multi-year ice concentration products

MYI concentration is an important indicator of abundance of older, thicker and rougher sea ice and is readily available from several resources including OSI SAF (Aaboe et al., 2017) and University of Bremen (Ye et al., 2016a). MYI can also be estimated from the NSIDC SIA product if we assume that bins with ice older than 1 year contain 100 % MYI and other bins contain 0 %. Total concentration of MYI from the SICCI product (with due consideration of SIC) was estimated as a sum of all multi-year fractions.

Intercomparison of MYI maps from four sources indicated that the general spatial distribution and temporal evolution is rather similar (Fig. 7): most of MYI is observed in the Canadian sector; a thin filament is exported through the Fram Strait; in some years the Beaufort Gyre advects MYI along the Beaufort Sea coast; and sleeves of enhanced MYI concentration are elongated towards the Laptev Sea.

Distinct features of MYI from NSIDC (Fig. 7a) include the following: the map is binary (ice/no ice); numerous discontinuities due to artifacts in the sea ice drift field are present; rather large lacunas filled mostly with FYI appear in the ice pack (e.g., 31 December 2015, Lincoln Sea); the MYI sleeves that extend towards the Laptev Sea are rather wide; and a wide gap between ice pack and the coastline is observed. The SICCI product (Fig. 7b) is continuous and smooth; expectedly it shows higher MYI concentrations in the central Arctic and near the Canadian coast. The OSI SAF MYI product (Fig. 7c) is also binary but the map is not discontinuous; generally it shows lower extent of MYI than other products, but in some periods the extent is much higher (e.g., on 29 March 2017 in the Beaufort Sea). The UB product is also continuous and apparently can resolve smaller spatial details but the atmospheric influence seems to reduce MYI concentration in the middle of the ice pack (e.g., along the 150^∘ meridian in 2013 or above the Yermak Plateau in 2015).

The time series of MYI areas are compared with regard to both seasonal dynamics (Fig. 8) and interannual variations (Fig. 9). Areas of ice age fractions of each age in the Arctic ocean are estimated individually from the NSIDC and SICCI products and plotted against time. For the NSIDC product the area of an ice age fraction is calculated as a number of all pixels corresponding to this ice age fraction multiplied by the area of the pixel. Similarly the MYI area from the OSI SAF product is estimated from the number of MYI pixels. For the SICCI product the area is estimated as the integral of the corresponding ice age fraction over all pixels. Similarly the MYI area from the UB product is estimated as the integral of MYI concentration over all pixels.

The seasonal variations of ice age fraction areas follow similar pattern for the NSIDC and SICCI products (Fig. 8): the sea ice area minimum (4 ×10⁶ km²) is observed in mid-September, and it is followed by a rapid growth of the FYI (blue in Fig. 8) until the beginning of the next year. A stable winter, when total area remain practically unchanged at level of 7–8 ×10⁶ km², is over in May and is followed by a rapid decrease of FYI. Unlike FYI, the area of MYI fractions cannot increase – it only gradually decreases due to ice convergence and melting. Maximum area of the 2YI is also observed in mid-September – at this point all FYI from the previous year is considered as the ice which has survived summer melting and it becomes 2YI. Other fractions of older ice increase their age by 1 year in the same fashion.

Figure 9 Interannual dynamics of area of sea ice age fractions (bars) and multi-year ice area (dots) derived from NSIDC product (left bars), from NERSC product (right bars), from UB product (black dots) and from OSI SAF (black stars).

The MYI area is shown in Figs. 8 and 9 as a demarcation of the blue and orange areas and as black and yellow dots. A comparison of MYI areas from the four products reveals several differences. The NSIDC MYI area decreases in two regimes: a relatively slow decrease during winter solely due to convergence and a faster decrease at the end of spring when melting reaches MYI. The SICCI MYI area decrease also accelerates towards spring but more gently – melting affects MYI already starting from the beginning of spring. The MYI areas from UB and OSI SAF correspond well to the SICCI algorithm but exhibit much more sporadic variability, including short periods of an inexplicable increase in MYI area. Although the MYI extent in the OSI SAF product is somewhat lower than of the other products (see maps in Fig. 7) the OSI SAF MYI area corresponds well to other estimates (Figs. 8 and 9) because concentration of MYI within the ice extent is assumed to be 100 %.

By the beginning of 2013 the SICCI product has only one fraction of MYI (2YI) because calculations were initiated in 2012 and all MYI was assigned to be 2YI. By the fourth annual cycle (in 2016) the SICCI product has five ice age fractions and the distribution can be better compared with the NSIDC product. Clearly, the fractions of the older ice in the SICCI product have lower area than in the NSIDC product.

The interannual variability can be seen on the plot of areas of the ice age fractions for 1 January of 2013, 2014, 2015 and 2016 in Fig. 9. The total ice area (entire height of the columns) is lower for the SICCI product because it is calculated as an integral ice concentration over all pixels whereas for the NSIDC it is calculated as sum of pixels within ice extent. The total MYI area is rather similar across the four products, but the area of older sea ice (≥2 years) fractions is lower by almost 20 % in the SICCI product.

4.3 Comparison of MYI with SAR

An independent validation of the SICCI MYI product was performed using a mosaic of synthetic aperture radar (SAR) images in HH polarization from Sentinel-1 A and B satellites acquired around 1 January 2016. MYI appears as brighter and more homogeneous texture on SAR images which allows one to draw an outline and compare it with outlines of MYI from the SICCI product as shown in Fig. 10. The SICCI MYI product corresponds very well to the SAR-derived MYI extent although it has much smoother boundaries due to low resolution of the input sea ice drift product (65 km). The low resolution is also the reason why SICCI MYI does not reach the coast, leaving a 50 km gap along the shoreline. The most considerable difference is observed in the Beaufort Sea where the SAR-based MYI is seen as a stretched “archipelago” of MYI “islands”. In this area the SICCI product exhibits erroneously smooth decline in MYI concentration and fills the gap between the MYI “archipelago” and the central Beaufort Sea.

Figure 10 Comparison of MYI extent derived by SAR visual interpretation (light blue) and from SICCI SIA product (red) on top of mosaic of Sentinel-1 SAR images for 1 January 2016. The thick red line shows 15 % threshold and thin red line shows 50 % threshold.

Figure 11 Comparison of sea ice age for 1 January 1985 initiated from ice parcels with various density (1 parcel per 12 × 12 km, 6 × 6 km, 3 × 3 km).

5 Discussion

The major advantage of the new algorithm is the ability to generate individual ice age fractions or, in other words, to provide a frequency distribution of the ice age in each pixel. This allows derivation of a single number characterizing SIA using one of the statistical methods presented above. Selection of the method depends on the preferred definition of the SIA and on the eventual application. For example, in our opinion, weighted average is the optimal method for graphical representation of SIA maps as it depicts the smooth nature of SIA distribution and accounts not only for the oldest ice fractions but also for the younger ones. However, for other applications, e.g., assimilation into models, a different method may be preferred. The SIA distribution can also be used to calculate MYI concentration for various applications, including calibration/validation of the PAMW-based MYI algorithms, improvement of the freeboard to sea ice thickness conversion for altimeter data and estimation of sea ice roughness for assimilation into models.

The motivation for implementing the Eulerian advection scheme in the new algorithm was to produce continuous and smooth spatial distributions of SIA fractions and also to prevent undersampling of the results. In the experiments with the NSIDC algorithm it was discovered that the density of starting points of the sea ice drift vectors indirectly influences the area of MYI. Too low density (e.g., 1 drifter in each 12 × 12 km box as designed at NSIDC) leads to undersampling of the ice age observations in the zone of strong divergence very early – only after a few years of propagation. With a two or four times higher number of initial points, undersampling is still experienced after 5 years of Lagrangian propagation. Observations that the Lagrangian approach introduces problems due to the dispersion of tracers were also previously reported in Rybak and Huybrechts (2003). The NSIDC method was run with three different initial densities of ice parcels: one drifter per 12 × 12 km or 6 × 6 km or 3 × 3 km box. The results show that lower initial densities result in a map of SIA with sparse presence of old ice fractions (Fig. 11a) and increasing the initial density of drift vectors leads to increase of concentration of the older ice (Fig. 11b and c). Consequently the area of older ice fractions ≥ 6 years) increases and the area of younger ice fractions (≤ 5 years) decreases when the initial density is increased. Impact of the increased initial density on increased MYI area is independent of the forcing and was demonstrated also for the OSI SAF ice drift field (not shown here).

The new sea ice drift product derived from AMSR2 also has a significant impact on the accuracy and reliability of the results. It accurately reproduces the sea ice circulation in the Arctic, is void of artifacts and produces homogeneous distribution of SIA fractions. It was discovered that the initial SID product is contaminated with high-frequency noise. Cleaning of SID with a median filter reduced small-scale deformations and prevented appearance of FYI discontinuities in artificial divergence zones. At the same time the low spatial resolution of SID (65 km) limits the algorithm from reproduction of fine-scale features in MYI distribution and may lead to over-smoothing in areas with high ice drift speed gradients (e.g., western Beaufort Sea).

The observed differences in MYI spatial distribution and total area are only marginal for the NSIDC and the SICCI products. This happens due to a short integration time: MYI spatial distribution is forced by observations on 15 September and then it is modified by an advection scheme (either Lagrangian or Eulerian) and melting (in case of SICCI product) during only 1 year. In contrast, the integration over longer times and accounting for the sea ice concentration have significant effect on the area of the older ice. Analysis of inter-annual dynamics (Fig. 9) shows that after 4 years of integration FYI and MYI areas in the NSIDC and SICCI products are almost the same (7.1 × 10⁶ and 6.9 × 10⁶ km²) but the area of the third-year and older ice in the SICCI product is almost 30 % smaller than in the NSIDC product. On average, the speed of annual multi-year ice decline (calculated as difference between MYI area in the beginning and in the end of the year) is 22 % higher for the SICCI product, which may also indicate that the actual residence time of sea ice is also shorter by one-fifth. It is challenging to separate the impact of the new sea ice drift product, scheme of advection and constraining by the SIC but more accurate estimates of long-term changes in the residence time will be possible when longer time series of the new SIA product become available.

The sea ice type algorithms based on PAMW satellite observations have higher spatial resolution than the ice-drift-based algorithms but inevitably suffer from the unaccounted atmospheric impact or melt ponds and have to be corrected using information on air temperature or ice motion (Ye et al., 2016a, b). The method in our paper has high potential for complementing the PAMW algorithms and production of high-resolution and consistent MYI estimations. Such combined procedure can be realized, for example, as a machine learning algorithm which is trained to output SICCI MYI using PAMW brightness temperatures (T_B) as input. The machine learning block can be realized as a simple polynomial regression or as a more complex neural network (Haykin, 1998). The training can be performed either on a per image basis, when only one mosaic of brightness temperatures is used as input and only one snapshot of MYI is used as output, or on a seasonal (or even satellite mission) basis, when information from many images is used for training. Alternatively, both SICCI MYI estimations and brightness temperatures can serve as input to classification algorithms, such as Bayesian classifier (Zhang, 2004) or support vector machines (Smola and Schölkopf, 2004).

6 Conclusions

We have developed a new algorithm for estimating sea ice age distribution using sea ice drift and concentration products. The algorithm is based on the Eulerian advection scheme which provides smooth distribution of the ice age parameters and prevents the undersampling problem that may occur when a Lagrangian tracking approach is used. Another advantage of the selected scheme is the ability to generate not just a single age characteristic but a distribution of sea ice age fractions. First, this allows for flexibility in choosing the ice age definition and application of a statistical measure to compute SIA and, second, this provides individual spatial distributions of ice age fractions that can be assimilated into models or used for ice type delineation. For example, concentration of multi-year ice can be computed as a sum of multi-year ice fractions and used for defining ice density and snow thickness for the ice thickness algorithms, ice roughness for the ice circulation models and so on.

The new algorithm is driven by the new sea ice drift products from OSI SAF, which is void of potential artifacts due to inclusion of autonomous ice drifter buoys. This leads to a more homogeneous distribution of ice age fractions over the Arctic Ocean. The algorithm is also constrained by the observed sea ice concentration from OSI SAF, which reduces fractions of old ice and, consequently, ice age by 20–30 %. It was applied to generate time series of daily sea ice age fraction product from October 2012 to October 2017. Comparisons with the NSIDC SIA time series indicate that the fractions of MYI in the new product melt faster during the year and after a spin-up time of 3 years the area of older ice in the SICCI product is almost 20 % lower than in the NSIDC product.