Temporal Means and Variability of Arctic Sea Ice Melt and Freeze Season Climate Indicators Using a Satellite Climate Data Record

Information on the timing of Arctic snow and ice melt onset, sea ice opening, retreat, advance, and closing, can be beneficial to a variety of stakeholders. Sea ice modelers can use information on the evolution of the ice cover through the rest of the summer to improve their seasonal sea ice forecasts. The length of the open water season (as derived from retreat/advance dates) is important for human activities and for wildlife. Long-term averages and variability of these dates as climate indicators are beneficial to business strategic planning and climate monitoring. In this study, basic characteristics of temporal means and variability of Arctic sea ice climate indicators derived from a satellite-based climate data record from March 1979 to February 2017 melt and freeze seasons are described. Our results show that, over the Arctic region, anomalies of snow and ice melt onset, ice opening and retreat dates are getting earlier in the year at a rate of more than 5 days per decade, while that of ice advance and closing dates are getting later at a rate of more than 5 days per decade. These significant trends resulted in significant upward trends for anomalies of inner and outer ice-free periods at a rate of nearly 12 days per decade. Small but significant downward trends of seasonal ice loss and gain period anomalies were also observed at a rate of −1.48 and −0.53 days per decade, respectively. Our analyses also demonstrated that the means of these indicators and their trends are sensitive to valid data masks and regional averaging methods.


Introduction
Since late 1978, passive microwave sensors have been providing a history of sea ice concentration from satellite data suitable for tracking climate change and variability in the polar regions.Sea ice extent (the area within the 15% concentration contour) and area (the area-integrated concentration) have long been considered as key climate indicators and have been included in numerous national and international climate assessment reports (e.g., [1,2]).However, these two parameters provide only limited information about the character of the sea ice; in addition, they have limited skill as indicators of future sea ice conditions, both seasonally and inter-annually [3].The day of opening (DOO) captures the start of the seasonal ice loss period, while the day of retreat (DOR) captures the end of the ice loss period [5].The 80% sea ice concentration (SIC) threshold typically marks the date when SIC declines toward the annual minimum [4,6] and is often used as the upper SIC threshold to define the marginal ice zone boundary.The 80% concentration threshold is also used by operational ice centers (e.g., U.S. National Ice Center [7]) as the boundary between the marginal ice zone and the pack ice.The difference between DOO and DOR is referred to as the seasonal loss of ice period (SLIP).Similarly, the day of advance (DOA) and day of closing (DOC) denote the beginning and ending of the seasonal ice gain period, respectively.The difference between the two is referred to as the seasonal gain of ice period (SGIP).The inner ice-free period (IIFP) is defined as the difference between DOA and DOR, which captures the period of open ocean when the grid cell ice fraction is less than 15%.The outer ice-free period (OIFP) is the difference between DOC and DOO, which captures both open ocean and ice transition periods.A schematic diagram is shown in Figure 1, which illustrates the temporal relationship between all dates (DOx) and periods.
Remote Sens. 2018, 10, x FOR PEER REVIEW 2 of 21 on the methods described in [4,5], using a long-term satellite-based climate data record.In this paper, we describe the temporal means and variability of the following indicators: snow and ice melt onset dates, sea ice opening, retreat, advance, and closing dates, seasonal ice loss and gain periods, inner and outer ice-free periods (as defined in Table 1 and described below.)The day of opening (DOO) captures the start of the seasonal ice loss period, while the day of retreat (DOR) captures the end of the ice loss period [5].The 80% sea ice concentration (SIC) threshold typically marks the date when SIC declines toward the annual minimum [4,6] and is often used as the upper SIC threshold to define the marginal ice zone boundary.The 80% concentration threshold is also used by operational ice centers (e.g., U.S. National Ice Center [7]) as the boundary between the marginal ice zone and the pack ice.The difference between DOO and DOR is referred to as the seasonal loss of ice period (SLIP).Similarly, the day of advance (DOA) and day of closing (DOC) denote the beginning and ending of the seasonal ice gain period, respectively.The difference between the two is referred to as the seasonal gain of ice period (SGIP).The inner ice-free period (IIFP) is defined as the difference between DOA and DOR, which captures the period of open ocean when the grid cell ice fraction is less than 15%.The outer ice-free period (OIFP) is the difference between DOC and DOO, which captures both open ocean and ice transition periods.A schematic diagram is shown in Figure 1, which illustrates the temporal relationship between all dates (DOx) and periods.Arctic snow and ice melt onset (MO) may be triggered by positive anomalies of air temperature and water vapor [8].The timing of MO on Arctic sea ice influences the amount of solar radiation absorbed by the ice-ocean system throughout the melt season by reducing surface albedos in the    Arctic snow and ice melt onset (MO) may be triggered by positive anomalies of air temperature and water vapor [8].The timing of MO on Arctic sea ice influences the amount of solar radiation absorbed by the ice-ocean system throughout the melt season by reducing surface albedos in the early spring [9].The timing of Arctic MO can be influenced by synoptic conditions [10].Significant trends of ~2-3 days/decade in early MO for the period of 1979-2011 were observed by [11].They pointed out that the variability in MO are largely driven by spring surface air temperature.Bliss et al. [9] compared two long-term time series of MO on Arctic sea ice derived from passive microwave brightness temperatures (Tbs) for the period of 1979-2012 and found large mean differences in trends between two different MO algorithms, with the largest uncertainty in thin-ice regions.Depending on the retrieval algorithms and correction methods, the Arctic MO decadal trends were found to range from −1.6 to −4.5 days/decade [9].
Stammerjohn et al. [12] identified trends in the DOR, DOA, and the duration of the sea ice season in the Kara and Barents Sea regions and for the East Siberian Sea, Chukchi Sea, and the western Beaufort Sea.Trends in the DOR and DOA for both regions ranged from 1 to 1.9 days per year (i.e., about 10-19 days per decade); with DOR occurring earlier in the year and DOA occurring later.The trends in DOR and DOA resulted in a shortening of the sea ice season by ~2.8 days per year (28 days per decade).Changes in the duration of the sea ice season were first reported by [13] who later identified a reduction in the length of the Arctic ice season by at least 5 days per decade [14].A strong relationship between the timing of ice retreat and that of maximum sea surface temperature (SST) was found by Steele and Dickinson [5].Stroeve et al. [15] investigated the predictability of DOA timing in response to timing of the prior DOR using several SIC thresholds and also found inverse correlations between DOR and DOA for most of the Arctic, as would be expected.However, some positive correlations were found in areas primarily in the eastern Arctic and along the ice edge near Franz Josef Land where the IIFP is short and can occur at varied times of the melt season from year to year.
The climatological means and temporal variability of these parameters will inform a variety of stakeholders.For example, the timing of melt onset and the date of sea ice retreat can provide useful information on the evolution of the ice cover through the rest of the summer, and thus be useful to efforts to improve seasonal sea ice forecasts (e.g., [15]).MO is also an indicator for summer minimum extent, since it plays a key role in the amount of solar energy absorption, and along with freeze-up, the overall energy balance of the Arctic climate system (e.g., [16,17]).The length of the open water season (as derived from retreat/advance dates) is important for human activities (navigability of waters, access to natural resources, hunting and transportation by indigenous populations) and for wildlife (e.g., polar bear access to food sources).The long-term time series will inform future planning of military, civilian, and commercial infrastructure (buildings, marine vessels).
Therefore, it is beneficial to provide users with the basic characteristics of these parameters from a consistent long-term time series.In this paper, we describe temporal mean and variability of Arctic climate indicators derived from a long-term satellite-based sea ice concentration climate data record.This is the first time that the basic statistics and temporal variability of all dates and periods are computed and described.By utilizing our integrated data set with all of these parameters, this study provides a new way to holistically look at the changing seasonality of the Arctic sea ice cover.
The paper is organized as follows.Section 2 outlines the datasets used in this study.Section 3 describes the temporal means and variability of these dates and periods.Discussions on sensitivity of different averaging choices and methods on the decadal trend of MO and on the means of dates and periods are provided in Section 4, followed by the conclusions in Section 5.

Datasets Outline
The daily passive microwave sea ice concentration climate data record (CDR) used in this study is created by a collaboration between NOAA (the National Oceanic and Atmospheric Administration) and NSIDC (National Snow and Ice Data Center).The SIC CDR leverages two well-established satellite products developed by NASA (the National Aeronautics and Space Administration) [18,19].It provides a single, authoritative sea ice concentration time series; all processing methods (including the software used to produce the data) are fully documented, reproducible, and publicly available [20].This is essential for establishing public confidence in the indicators.
The sea ice opening, retreat, advance, and closing dates are derived from the merged Goddard SIC fields from the daily CDR data files for each ice year, starting from 1 March each year and ending on February 28 the following year.The SIC fields are on the NSIDC polar stereographic grid with nominal 25 km × 25 km grid cells (see [21] for a detailed description of variables in the CDR data files).Extending the ice year to the end of February of the following year captures dates of ice closing that may occur beyond 31 December.
SIC is smoothed with a 7-day boxcar running mean prior to computing dates to help filter out potentially spurious ice variability associated with synoptic events.Prior to July 1987, the SIC data are only available every other day and a linear interpolation was used to fill the gaps.For consistency over a melt season, DOA and DOC are only derived for grid cells where a corresponding DOR and DOO are identified.In some cases, a valid DOA or DOC may not be observed before the end of our defined ice year.This typically occurs along the southern ice edge where ice extent has not yet reached its maximin on 28 February.
Snow and ice MO dates are calculated directly from brightness temperatures, the same source data used for the SIC, which provides consistency across the indicators.The advanced horizontal range algorithm (AHRA) MO dates give the initial melting date of the snow and ice after which freeze/thaw cycling may continue to occur [22,23].

Computing Regional Means
Sea ice retreat has been extreme in recent years, but this was not the case in the earlier years of the satellite era.This creates a quandary for generating time series of regional mean parameters with respect to the area over which to average.There are two choices: (i) average in each year only over the area for which a parameter is valid for ALL years (i.e., the minimal intersection of all years), referred to as Case core or (ii) average in each year over the full area valid for that parameter in that year, referred to as Case allv.
An analysis is carried out to examine both strategies on their impact to trends using the [23] MO data.The basic statistics are captured in Table 2.While the cell number remains the same at 17,058 for Case core, by definition the total valid date cells for a given year (Case allv) range from 20,883 to 24,717.The difference in mean grid cell numbers between these two cases is 6221, which is about one-third of the valid date cell number for Case core.There is a small but significant downward trend in the valid date cell numbers for Case allv (top panel in Figure 2).This is driven by the negative trend in March sea ice extents [19].The MO dates are calculated for sea ice grid cells that approximate the annual maximum sea ice extent.The start of melt season sea ice extent is set for the first two days with SIC data from the CDR for March (note that this means March 1st and 3rd or March 2nd and 4th for SMMR years).If the SIC is 50% or greater for one or both of the first two days in March, a MO date is calculated for that cell.Since the extent of the sea ice edge changes from year to year, it will result in a different number of grid cells where a MO date was calculated each year.As a result of the downward trend in the annual maximum sea ice extent, the numbers of grid cells with SIC 50% or greater for one or both of the first two days in March has been reduced.

Results
The spatial distributions of long-term averages of dates and periods, their standard deviation (STD), and significant trends are first described.Then, the decadal trends of these dates and periods averaged over the whole Arctic region are computed and examined.Unless otherwise mentioned, significant trends hereinafter refer to those that are significant at the 95% confidence level.The pole hole diameter also changes depending on the SIC CDR pole hole mask which uses the SMMR, SSM/I, or SSMIS pole holes (Table 3).There are two signals represented in the Case allv cell counts: (1) the decrease in March SIE, and (2) a step change between 1987 and 1988 when the source data shifts from SMMR to SSM/I and the pole hole diameter changes; again, between 2007 and 2008 from SSM/I to SSMIS (Table 3).This downward trend in Case allv would be more extreme if the size of the pole hole was constant over the entire time series.The trends in the middle panel of Figure 2 are from the time series of regional averages of the valid melt dates within that region for both allv and core cases (hereinafter referred to as Method I).The mean that goes into the calculation in each case contains both spatial and temporal averages.In this approach, the mean remains the same in both temporal and space domains.Noticeable differences in both time mean and decadal trends can be observed with a mean bias of −11.2 days and 3.08 days/decade, respectively (Figure 2 and Table 4).The mean bias is largely due to the removal of the early melting sea ice area along the shifting edge of the sea ice pack, which are filtered from the Case core regional averages resulting in the later MO bias relative to Case allv.The trends in the bottom panel of Figure 2, on the other hand, are from the time series of regional means of MO departures from its own long-term mean (1979-2016) at each grid cell (hereinafter referred to as Method II).As the long-term temporal means are spatially varied, this approach will remove the impact of the spatially varied time mean and isolate the temporal variations.The long-term mean of the MO anomalies for both cases are very close to each other, which are at or close to 0 (Table 4).There is also only a slight difference in the linear decadal trends in both cases-the mean bias is only 0.7 day/decade (Table 4), which represents the effect by the transient valid date cells.
The approach of utilizing Method II and Case allv is adopted in the rest of this paper.Unless otherwise mentioned, all the linear decadal trends are computed from the time series of regional averages of anomalies (namely, departures from their long-term temporal averages) at individual grid cells.

Results
The spatial distributions of long-term averages of dates and periods, their standard deviation (STD), and significant trends are first described.Then, the decadal trends of these dates and periods averaged over the whole Arctic region are computed and examined.Unless otherwise mentioned, significant trends hereinafter refer to those that are significant at the 95% confidence level.3a).The latest melting region of the Arctic is in the North American sector of the central Arctic to the north of the Canadian Arctic Archipelago and Greenland where mean surface air temperatures are generally coldest and the ice is predominantly thicker multiyear ice.MO in this region begins in early to mid-June.MO in the eastern Arctic is typical during the month of May.Interannual variability in MO is generally high with standard deviations of the date of MO as large as 40 days in parts of the North Atlantic sector (Figure 3b).For much of the Arctic, STD in the dates of MO exceed one month as a result of variable spring weather conditions which initiate MO.

Spatial Distribution of Mean and STD of MO
The mean date of MO for the Arctic during the 1979-2016 period is latitudinally dependent on earlier MO occurring in the southernmost latitudes and along the peripheral sea ice edge during March and April (Figure 3a).The latest melting region of the Arctic is in the North American sector of the central Arctic to the north of the Canadian Arctic Archipelago and Greenland where mean surface air temperatures are generally coldest and the ice is predominantly thicker multiyear ice.MO in this region begins in early to mid-June.MO in the eastern Arctic is typical during the month of May.Interannual variability in MO is generally high with standard deviations of the date of MO as large as 40 days in parts of the North Atlantic sector (Figure 3b).For much of the Arctic, STD in the dates of MO exceed one month as a result of variable spring weather conditions which initiate MO.
The spatial patterns in the MO date mean and STD compare similarly with those from other studies of pan-Arctic MO from passive microwave observations for most of the Arctic [6,11,24] with the exception that the AHRA MO dates occur earlier in the year than MO dates from these other studies.Bliss et al. [9] investigated the differences between the AHRA MO and the earliest MO dates from the Markus et al. [6] method, finding a mean difference of ~20 days due to differences in the sensitivity of the frequency channels used and algorithm methodology.A threshold method in the Markus et al. [6] algorithm produces MO dates that are more similar to the DOO for some locations in the peripheral ice regions; thus, we used the AHRA MO dates which are earlier and more closely tied to the spring atmospheric conditions that initiate melt as used by others [8,10].Regionally, other methods have been used to identify MO dates in the Canadian Arctic using synthetic aperture radar (SAR) imagery (e.g., [25]).Passive microwave MO dates from both the AHRA method (Figure 3) and the Markus et al. [6] method were on average 7 days earlier than the SAR MO dates in this region [9].The later SAR MO dates are likely related to the local early morning image acquisition time when surface air temperatures are cooler than the daily mean air temperatures, and the less frequent repeat observation time (2-4 days) of the SAR imagery.

Spatial Distribution of Mean and STD of DOx and Periods
DOO and DOR are generally earlier in the south and later in the north (Figure 3a).An interesting exception is in Hudson Bay, which melts out from north-to-south due to prevailing northwesterly winds that push ice south from northern Hudson Bay in late spring/early summer [26].Another exceptional area is the southern Nansen Basin (to the north of Svalbard and eastward from there), The spatial patterns in the MO date mean and STD compare similarly with those from other studies of pan-Arctic MO from passive microwave observations for most of the Arctic [6,11,24] with the exception that the AHRA MO dates occur earlier in the year than MO dates from these other studies.Bliss et al. [9] investigated the differences between the AHRA MO and the earliest MO dates from the Markus et al. [6] method, finding a mean difference of ~20 days due to differences in the sensitivity of the frequency channels used and algorithm methodology.A threshold method in the Markus et al. [6] algorithm produces MO dates that are more similar to the DOO for some locations in the peripheral ice regions; thus, we used the AHRA MO dates which are earlier and more closely tied to the spring atmospheric conditions that initiate melt as used by others [8,10].Regionally, other methods have been used to identify MO dates in the Canadian Arctic using synthetic aperture radar (SAR) imagery (e.g., [25]).Passive microwave MO dates from both the AHRA method (Figure 3) and the Markus et al. [6] method were on average 7 days earlier than the SAR MO dates in this region [9].
The later SAR MO dates are likely related to the local early morning image acquisition time when surface air temperatures are cooler than the daily mean air temperatures, and the less frequent repeat observation time (2-4 days) of the SAR imagery.
3.1.2.Spatial Distribution of Mean and STD of DOx and Periods DOO and DOR are generally earlier in the south and later in the north (Figure 3a).An interesting exception is in Hudson Bay, which melts out from north-to-south due to prevailing northwesterly winds that push ice south from northern Hudson Bay in late spring/early summer [26].Another exceptional area is the southern Nansen Basin (to the north of Svalbard and eastward from there), which is strongly controlled by warm water advection from the West Spitzbergen Current [27].It takes a long time for this current to cool in the late summer and fall, resulting in much-delayed DOO and DOR.
The strongest gradient of early-to-late dates is found farther south in the Pacific Sector (Okhotsk-Bering-Chukchi-Beaufort), while the mean ice edge is more northerly in the Atlantic Sector (Barents-Greenland-Baffin-Kara).Within the Arctic Ocean, particularly early DOO and DOR occur in the eastern Beaufort Sea, Chukchi Sea, and Laptev Sea (Figure 4a), owing mostly to late spring offshore wind forcing [4].
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 21 which is strongly controlled by warm water advection from the West Spitzbergen Current [27].It takes a long time for this current to cool in the late summer and fall, resulting in much-delayed DOO and DOR.
The strongest gradient of early-to-late dates is found farther south in the Pacific Sector (Okhotsk-Bering-Chukchi-Beaufort), while the mean ice edge is more northerly in the Atlantic Sector (Barents-Greenland-Baffin-Kara).Within the Arctic Ocean, particularly early DOO and DOR occur in the eastern Beaufort Sea, Chukchi Sea, and Laptev Sea (Figure 4a), owing mostly to late spring offshore wind forcing [4].
Even though DOR is relatively early in the Laptev Sea, SSTs do not warm as much as in the Beaufort Sea with comparable early DOR, owing to the Laptev's more northerly position [5].As a result, DOA and DOC are not particularly late in this area, in contrast to the anomalously late freezeup seen in the Beaufort Sea (Figure 4a).The Chukchi has a particularly late freeze-up, owing in part to warm southerly ocean currents that move into this area during summer (i.e., in addition to the effect of solar radiation absorption into open water) [28].
(a) The mean SLIP has strong spatial variability compared to the mean SGIP (Figure 5b).In the south, a relatively short SLIP is likely forced mostly by strong downward atmospheric heat fluxes.Just to the north, however, wind forcing blows the ice pack around while the atmospheric heat fluxes are still relatively weak.This leads to a highly variable, "noisy" ice concentration time series that allows for a long SLIP (e.g., [4]).Finally, in the far north, advection of warm ocean water (which very efficiently and quickly melts ice) plays a dominant role in sea ice retreat [29].
In contrast, the mean SGIP has relatively small spatial variability (Figure 5b).The period of ice gain is more thermodynamically controlled relative to the SLIP, i.e., vast areas of ice can freeze up all at once (even if the ice is thin).However, there can be "noisy" ice concentration time series during this period as well, as surface cooling and new ice growth is followed by mixing up of subsurface water warmed during the previous summer [30].This is an area of active research.
High DOO and DOR variance is found in the same area with long SLIP, likely owing to windforced variability (Figure 5a).High variance in all DOx fields is found in the southern Nansen Basin, Even though DOR is relatively early in the Laptev Sea, SSTs do not warm as much as in the Beaufort Sea with comparable early DOR, owing to the Laptev's more northerly position [5].As a result, DOA and DOC are not particularly late in this area, in contrast to the anomalously late freeze-up seen in the Beaufort Sea (Figure 4a).The Chukchi has a particularly late freeze-up, owing in part to warm southerly ocean currents that move into this area during summer (i.e., in addition to the effect of solar radiation absorption into open water) [28].
The mean SLIP has strong spatial variability compared to the mean SGIP (Figure 5b).In the south, a relatively short SLIP is likely forced mostly by strong downward atmospheric heat fluxes.Just to the north, however, wind forcing blows the ice pack around while the atmospheric heat fluxes are still relatively weak.This leads to a highly variable, "noisy" ice concentration time series that allows for a long SLIP (e.g., [4]).Finally, in the far north, advection of warm ocean water (which very efficiently and quickly melts ice) plays a dominant role in sea ice retreat [29].

Spatial Distribution of Decadal Trends of MO
Local trends in the date of MO over the 1979-2016 period are computed from the time series of MO dates at each grid cell (Figure 3c).Only statistically significant trends are shown.The majority of the Arctic is dominated by negative trends in MO date indicating a shift towards MO earlier in the year.
Statistically significant positive trends in MO date are present in a few locations in the Bering Sea where a positive trend in winter SIE was observed by [31].Cooler temperatures and an increase in northwesterly winds during the 2003-2011 period have been attributed to the extensive sea ice cover in this region [32] which would contribute to the observed later MO date trends.Other small areas with positive trends occur largely along the coastlines in northern Hudson Bay, the Canadian Arctic Archipelago, and the Sea of Okhotsk where land contamination in the passive microwave data [33], likely causing spurious MO dates [34].In contrast, the mean SGIP has relatively small spatial variability (Figure 5b).The period of ice gain is more thermodynamically controlled relative to the SLIP, i.e., vast areas of ice can freeze up all at once (even if the ice is thin).However, there can be "noisy" ice concentration time series during this period as well, as surface cooling and new ice growth is followed by mixing up of subsurface water warmed during the previous summer [30].This is an area of active research.
High DOO and DOR variance is found in the same area with long SLIP, likely owing to wind-forced variability (Figure 5a).High variance in all DOx fields is found in the southern Nansen Basin, likely owing to interannual variability in the heat transport of Atlantic Waters in the West Spitzbergen Current.

Spatial Distribution of Decadal Trends of MO
Local trends in the date of MO over the 1979-2016 period are computed from the time series of MO dates at each grid cell (Figure 3c).Only statistically significant trends are shown.The majority of the Arctic is dominated by negative trends in MO date indicating a shift towards MO earlier in the year.
Statistically significant positive trends in MO date are present in a few locations in the Bering Sea where a positive trend in winter SIE was observed by [31].Cooler temperatures and an increase in northwesterly winds during the 2003-2011 period have been attributed to the extensive sea ice cover in this region [32] which would contribute to the observed later MO date trends.Other small areas with positive trends occur largely along the coastlines in northern Hudson Bay, the Canadian Arctic Archipelago, and the Sea of Okhotsk where land contamination in the passive microwave data [33], likely causing spurious MO dates [34].

Spatial Distribution of Significant Trends of DOx and Periods
Ice opening and loss dates (DOO and DOR) are generally trending earlier, while ice gain and closing dates (DOA and DOC) are trending later, with all areal-mean trends roughly equal in absolute magnitude at ~7-10 days/decade (Figure 6a).This analysis only computes trends at locations with a DOx in at least 10 years over the entire record, so (for example) this does not show the very large changes occurring in the far north in the past few years.Trends in DOx variables are generally weak in subpolar seas and strongest in the marginal seas of the Arctic Ocean.The eastern Beaufort Sea has relatively weak trends in all DOx variables over this long-term time series, although more recent years show stronger trends [4,32].
Trends in the SLIP and SGIP are generally small, although some significant and spatially variable SLIP trends are found in the western Arctic Ocean marginal seas (Figure 6b).Strong positive trends are evident in IIFP and OIFP.Therefore, the Arctic Ocean is becoming less icy over the year.Ice opening and loss dates (DOO and DOR) are generally trending earlier, while ice gain and closing dates (DOA and DOC) are trending later, with all areal-mean trends roughly equal in absolute magnitude at ~ 7-10 days/decade (Figure 6a).This analysis only computes trends at locations with a DOx in at least 10 years over the entire record, so (for example) this does not show the very large changes occurring in the far north in the past few years.Trends in DOx variables are generally weak in subpolar seas and strongest in the marginal seas of the Arctic Ocean.The eastern Beaufort Sea has relatively weak trends in all DOx variables over this long-term time series, although more recent years show stronger trends [4,32].
Trends in the SLIP and SGIP are generally small, although some significant and spatially variable SLIP trends are found in the western Arctic Ocean marginal seas (Figure 6b).Strong positive trends are evident in IIFP and OIFP.Therefore, the Arctic Ocean is becoming less icy over the year.The decadal trend of Arctic MO departures is −5.23 days/decade, which is significant at the 99% confidence level (Table 4).As the snow and ice melt onset dates move earlier from its long-term average state at the rate of more than 5 days per decade, it will potentially influence the Arctic iceocean system by influencing the amount of solar radiation.

Arctic DOx and Periods Decadal Trends
The mean values of Arctic DOx and periods are listed in Table 5 for four different valid date mask cases.In addition to Case allv and Case core, two other cases are considered here.Case mask includes only cells with valid DOR.On the other hand, Case comm requires dates valid for all DOx

Decadal Arctic MO Trends
The decadal trend of Arctic MO departures is −5.23 days/decade, which is significant at the 99% confidence level (Table 4).As the snow and ice melt onset dates move earlier from its long-term average state at the rate of more than 5 days per decade, it will potentially influence the Arctic ice-ocean system by influencing the amount of solar radiation.

Arctic DOx and Periods Decadal Trends
The mean values of Arctic DOx and periods are listed in Table 5 for four different valid date mask cases.In addition to Case allv and Case core, two other cases are considered here.Case mask includes only cells with valid DOR.On the other hand, Case comm requires dates valid for all DOx for the same ice year, and therefore contains only those with a complete ice opening, retreat, advance, and closing cycle.The means from Case core include interannual variability of cells that are valid in all years but do not include the contribution by the cells that are transient over the entire record period.Case comm captures the interannual variability of complete ice cycles.The fact that the results computed for Case core are very different from that for Case allv and Case comm implies that the transition cells and interannual variability of ice cycles are changing the basic characteristics of DOx, and therefore periods.Given the rapidly changing sea ice coverage in the Arctic, regularly monitoring these climate indicators and frequently updating their basic characteristics are necessary.
Sea ice opening and closing dates (DOO and DOC, respectively) are most sensitive to the valid date masks, and result in the OIFP (Table 5).While the results are fairly similar between Case mask and Case comm, the DOO mean for Case allv is about 10 days more and the DOC mean is 14 days less than that from Case mask or Case comm.In Case allv, it is especially counter-intuitive for the DOC mean to be less than the DOA mean but it is actually plausible.Additional discussion is included in Section 4.
All DOx anomalies exhibit decadal trends that are significant at the 99% confidence level.Both DOO and DOR anomalies show significant downward trends of about −5.29 and −6.27 days/decade, respectively (Figure 7 and Table 6).These downward trends imply that the opening and retreat dates are getting earlier in the year.The advance and closing dates, on the other hand, all display significant upward trends, ranging from 5.63 to 6.52 days/decade (Figure 7 and Table 6).These upward trends imply that the advance and closing dates are moving back in the year.* indicates statistical significance at the 99% confidence level.The combination of these earlier opening and retreat dates and later advance and closing dates leads to significant upward trends for both inner and outer ice-free period anomalies with decadal trend values of nearly 12 days/decade (Figure 7 and Table 6).The seasonal loss and gain of ice period anomalies display small downward significant trends of about -1.48 and -0.53 days/decade, respectively (Figure 7 and Table 6).All trends are significant at the 99% confidence level.
The seasonal gain of ice period, namely, SGIP, tends to be shorter than that of the seasonal loss of ice period, namely, SLIP (see Table 5).This implies that the melting process is much slower than that of the freeze-up process.
The total numbers of grid cells with valid dates for the whole Arctic region experienced significant upward linear decadal trends, ranging from 4.32%/decade for DOO to 7.41%/decade for DOA (Figure 8).The percentage increase is fairly similar for the ice opening and closing, namely, DOO and DOC, while the same is true for ice retreat and advance, namely, DOR and DOA.Therefore, the results show a bigger, increased trend in the total numbers of grid cells that are undergoing ice retreat and advance over the whole Arctic.On a year-to-year basis, the change is more variable and the increase is less distinct for all dates.
Remote Sens. 2018, 10, x FOR PEER REVIEW 16 of 21 The total numbers of grid cells with valid dates for the whole Arctic region experienced significant upward linear decadal trends, ranging from 4.32%/decade for DOO to 7.41%/decade for DOA (Figure 8).The percentage increase is fairly similar for the ice opening and closing, namely, DOO and DOC, while the same is true for ice retreat and advance, namely, DOR and DOA.Therefore, the results show a bigger, increased trend in the total numbers of grid cells that are undergoing ice retreat and advance over the whole Arctic.On a year-to-year basis, the change is more variable and the increase is less distinct for all dates.The mean seasonal evolution of Arctic sea ice cover can be summarized by relating the variability in the sea ice indicators during the spring through autumn summer sea ice loss cycle from year to year (Figure 9).Interannual variability of all dates are noticeable but do not seem to be highly correlated with significant years in sea ice history, such as 1998, 2007, and 2012.Additional in-depth study will be beneficial in shedding more light on this particular issue.The mean seasonal evolution of Arctic sea ice cover can be summarized by relating the variability in the sea ice indicators during the spring through autumn summer sea ice loss cycle from year to year (Figure 9).Interannual variability of all dates are noticeable but do not seem to be highly correlated with significant years in sea ice history, such as 1998, 2007, and 2012.Additional in-depth study will be beneficial in shedding more light on this particular issue.

Discussion
The total number of valid dates cells varies on a year-to-year basis.As mentioned previously, there are two choices for computing regional averages: (i) average in each year only over the area for which a parameter is valid for ALL years (i.e., the minimal intersection of all years), or (ii) average in each year over the full area valid for that parameter in that year.Large difference of more than 11 days was found in regional averages of MO dates between these two different averaging choices.So was their decadal trend difference of about 3 days per decade, which is as high as the trend computed from the time series obtained from the first averaging choice (see Figure 2 and Table 4).However, the differences are greatly reduced if the long-term means are removed at each valid date cell before the regional averaging is carried out.The biases are now only 0.09 day between these time series of MO

Discussion
The total number of valid dates cells varies on a year-to-year basis.As mentioned previously, there are two choices for computing regional averages: (i) average in each year only over the area for which a parameter is valid for ALL years (i.e., the minimal intersection of all years), or (ii) average in each year over the full area valid for that parameter in that year.Large difference of more than 11 days was found in regional averages of MO dates between these two different averaging choices.So was their decadal trend difference of about 3 days per decade, which is as high as the trend computed from the time series obtained from the first averaging choice (see Figure 2 and Table 4).However, the differences are greatly reduced if the long-term means are removed at each valid date cell before the regional averaging is carried out.The biases are now only 0.09 day between these time series of MO anomalies and 0.7 day per decade for their decadal trends.While the trend computed for the first averaging choice remains the same for both time series of MO and MO anomaly, the trend computed from the second averaging choice using the MO anomalies is very close to that from the first averaging choice.This demonstrates that the spatial-varied mean state of valid dates has a big impact on the trend analysis.If the second averaging choice is preferred to focus on temporal variability, one should use the anomaly fields to compute regional time series.This may, in part, explain the differences in trend values of dates and periods between our results and the previous results, such as those from [12,14].
After the initial melting date, additional freeze/thaw cycles may occur.Then, the question is whether continuous melt onset (CMO) will be a better indicator.The passive microwave CMO data available rely on a sea ice concentration threshold in cases when a clear melting signal cannot be determined from the brightness temperatures [6].In these cases, the algorithm for continuous melt identifies the day when the SIC falls below 80% for the last time before the location remains ice-free until autumn freeze-up.The threshold is primarily used to identify continuous melt onset in the peripheral ice regions [9] where we obtain our dates and periods of retreat/advance.This threshold conflicts with the DOO indicator, which also uses an 80% SIC threshold.In this case, we chose to use the initial MO date to identify the earliest changes in the sea ice cover at the beginning of the melt season.The earlier MO dates are more closely tied to the spring atmospheric conditions that initiate melt (as used by Mortin et al. [8] and Liu and Schweiger [10]) and therefore give additional information about the sea ice melt and retreat period that the continuous MO date would not provide.Bliss et al. [9] have compared the initial, early and continuous MO dates.The annual evolutions of those MO dates for our record periods are included in Figure 9 as reference.
Because the periods are inter-related with sea ice dates, calculation of a period mean includes only date cells that are valid to both the dates used to compute the period.Due to the fact that not all valid date cells undergo a complete melt, retreat, advance, and freeze-up cycle, it is expected to see some sensitivity of the resultant mean values to which valid date cells are used to compute the means.For Case allv, where all valid date cells are included when calculating the Arctic means of individual dates and periods, the DOC mean is 309.3, which is actually short of the DOA mean of 309.8, while the SGIP mean is 14.6 days (Table 5).This may seem to be counter-intuitive but it is in fact plausible as it is not unusual to have large areas of very early DOC in the north, namely, right after the sea ice minimum.This happens because large areas of the interior drop below 80% at the end of summer (namely, very late DOO), and then close right back up again just after the sea ice minimum (namely, very early DOC).The calculation of SGIP mean uses only those date cells that are valid for both DOA and DOC, which will mask out those cells.Additional study is needed to confirm this, which is beyond the scope of this study.We have, however, considered two other cases to help demonstrate the potential validity of our speculation.Case mask is when the cells are masked out by that of DOR.Case mask implies that only cells that experienced ice retreat are selected.Case comm is only used for cells, of which all dates are valid for the year.Case comm implies that all the selected cells have undergone the complete melt, retreat, advance, and freeze-up cycle.The total records for computing DOC means are reduced in both Case mask and Case comm, compared to that of Case allv (Table 5).The DOC means for Case mask and Case comm are 323.1, which is 13.3 and 14.6 days more than the DOA means, respectively.

Conclusions
Temporal means and variability of a set of climate indicators that track the seasonal evolution of the Arctic snow and ice cover have been described.These climate indicators, derived from a well-managed and documented satellite CDR, include snow and ice melt onset dates, sea ice opening and retreat dates, sea ice advance and closing dates, seasonal ice loss or gain periods, inner and outer ice-free periods.
Significant trends calculated over the period of March 1979-February 2017 are observed for all dates (MO, DOO, DOR, DOA, and DOC).The dates of melt onset, opening, and retreat are getting earlier in the year at a rate of more than 5 days per decade.The dates of ice advance and closing are getting later in the year at a rate of more than 5 days per decade.As the result, both inner ice-free period (i.e., potential open ocean) and outer ice-free period show a significant upward trend of nearly 12 days per decade.Seasonal ice loss period (SLIP) and seasonal ice gain period (SGIP) exhibit a slight but significant downward trend of −1.48 days per decade and −0.53 days per decade, respectively.All trends of dates and periods are significant at the 99% confidence level.
If these significant trends persistent, there will be potential impact on human activities and for wildlife due to the increased length of ice-free periods.Earlier melt onset and prolonged open water reduce surface albedos which could have a large influence on the total amount of solar radiation absorbed by the Arctic ice-ocean system [24].
Our analysis has demonstrated the sensitivity of computing long-term means of dates and periods, and that of the MO trends to regional averaging methods and valid cells masks due to the transient valid date cells.This is exacerbated by recent accelerated Arctic sea ice melting.It is our hope to start a community-wide dialogue to standardize how they should be computed for the purpose of climate monitoring of these indicators.
Sea ice coverage varies on the regional scales (e.g., [35]).Variability in trends for dates and periods does occur on the regional scales (see Figure 6; also, [6,12,36]).The regional variability of the climate indicators will be examined in more detail by [37], using the same datasets utilized in this paper.
The retrieval algorithms of passive microwave sea ice concentrations are sensitive to emissivity and surface temperature of sea ice, atmospheric effects, melt ponds, and thin ice, etc.They are also prone to land contamination in the coastal areas.Large uncertainty in SIC retrievals are found in the marginal ice zone, especially in summer, due to atmospheric and wind roughness efforts of open water areas [38].(See [39] for a comprehensive review of accuracy and precision of sea ice concentration algorithms in various sea ice and atmospheric conditions and inter-comparison of the SIC retrieval algorithms.)Methods have been developed to minimize the impact of those error sources [20].The error sources and uncertainties associated with SIC algorithms will nonetheless adversely affect the accuracy and precision of the DOx algorithms.
We here use a consistent and authoritative SIC data set (the CDR), which ensures the best possible inputs to our DOx algorithms.The existing DOx algorithms differ slightly among the small (but growing) number of researchers investigating the seasonal timing of sea ice retreat and advance.Going forward, an inter-comparison of these algorithms, including an in-depth study of a few significant years in sea ice history, for example, 1998, 2007, 2012, etc., would be useful.
Author Contributions: The Arctic sea ice spring retreat and fall freeze-up dates products are developed by M.S. and S.D. S.D. carried out the analysis and produced Figures 4-6. A. Bliss carried out the analysis and produced Figures 3 and 9. G.P. carried out the rest of analyses and produced all other figures and all tables.All authors have contributed to or reviewed the draft and provided beneficial feedback.
Funding: This work was funded by NASA's Climate Assessment Products and Indicators Program under the grant NNX16AK43G.Acknowledgments: J.A. Miller provided the PMW algorithm EMO and CMO time series plotted in Figure 9. Andrew Ballinger, Tom Maycock, and Jason Yu contributed to the design of Figure 1.Jessica Griffin has reviewed Figure 1 with beneficial comment.She also polished the graphical abstract.We thank Russell Vose, Matthew Menne, and David Babb for beneficial discussions.Comments from the Remote Sensing anonymous reviewers have improved the clarity and readability of the paper.

Conflicts of Interest:
The authors declare no conflict of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Figure 1 .
Figure 1.A schematic diagram shows the sea ice concentration (SIC) threshold values for DOO, DOR, DOA and DOC.It also displays the relationships between dates and periods, namely, SLIP, SGIP, IIFP, and OIFP, for a sea ice opening, retreat, advance, and closing cycle.It is for demonstration only and not to scale.See Table1for the acronyms and their definitions.

Figure 1 .
Figure 1.A schematic diagram shows the sea ice concentration (SIC) threshold values for DOO, DOR, DOA and DOC.It also displays the relationships between dates and periods, namely, SLIP, SGIP, IIFP, and OIFP, for a sea ice opening, retreat, advance, and closing cycle.It is for demonstration only and not to scale.See Table1for the acronyms and their definitions.

Figure 2 .
Figure 2. The number of valid snow melt onset date grid cells (×1000) (top), snow melt onset date (middle), and the departure (bottom) for the whole Arctic.Pink denotes Case allv where all the valid date cells are included in the averaging.Green denotes Case core where only the cells have a valid date for the whole record period of (1979-2016).Each ice year runs from the beginning of March to the end of February of the following year.The solid lines with circles are the time series of the variable.The thin dotted lines are for the time means.The thick dashed lines denote the decadal trend lines.

Figure 2 .
Figure 2. The number of valid snow melt onset date grid cells (×1000) (top), snow melt onset date (middle), and the departure (bottom) for the whole Arctic.Pink denotes Case allv where all the valid date cells are included in the averaging.Green denotes Case core where only the cells have a valid date for the whole record period of (1979-2016).Each ice year runs from the beginning of March to the end of February of the following year.The solid lines with circles are the time series of the variable.The thin dotted lines are for the time means.The thick dashed lines denote the decadal trend lines.

3. 1 .
Spatial Distributions of Climatological Means and STDs of Dates and Periods 3.1.1.Spatial Distribution of Mean and STD of MO The mean date of MO for the Arctic during the 1979-2016 period is latitudinally dependent on earlier MO occurring in the southernmost latitudes and along the peripheral sea ice edge during March and April (Figure

Figure 3 .
Figure 3.The mean (a), standard deviation (STD) (b), and linear decadal trend of snow and melt onset dates for the period of 1979-2016 (c).

Figure 3 .
Figure 3.The mean (a), standard deviation (STD) (b), and linear decadal trend of snow and melt onset dates for the period of 1979-2016 (c).

Figure 5 .
Figure 5. (a) Same as Figure 4a except for standard deviation (STD).The time series at each grid cell has been de-trended and the minimum number of the records for STD calculation is 10 years; (b) Same as (a) except for SLIP, IIFP, OIFP, and SGIP.

Figure 6 .
Figure 6.(a) Same as Figure 4a except for significant decadal trends (days/decade).Only those that are significant at the 95% confidence level are plotted.The minimum number of records for the trend calculation is 10 years; (b) Same as (a) except for SLIP, IIFP, OIFP, and SGIP.

Figure 7 .
Figure 7. Time series of DOO, DOR, DOA, DOC, SLIP, SGIP, IIFP, and OIFP anomalous (days, red circles with solid black line) and their linear regression trend lines (thick green dashed lines).The dash-dotted lines are one standard deviation of each time series.The decadal trends in red are significant at the 99% confidence level using Student's t-test.MN and STD denote the mean and standard derivation of the time series.

Figure 7 .
Figure 7. Time series of DOO, DOR, DOA, DOC, SLIP, SGIP, IIFP, and OIFP anomalous (days, red circles with solid black line) and their linear regression trend lines (thick green dashed lines).The dash-dotted lines are one standard deviation of each time series.The decadal trends in red are significant at the 99% confidence level using Student's t-test.MN and STD denote the mean and standard derivation of the time series.

Figure 8 .
Figure 8.The total number of grid cells (×1000) with valid DOx dates, namely DOO, DOR, DOA, and DOC for the whole Arctic region.

Figure 8 .
Figure 8.The total number of grid cells (×1000) with valid DOx dates, namely DOO, DOR, DOA, and DOC for the whole Arctic region.

Figure 9 .
Figure 9. Mean annual evolution of the sea ice melt season for the Northern Hemisphere (March 1979-February 2017) for Case mask.Colors for the shaded bars define the mean SLIP, IIFP, SGIP periods.The full span of the shaded bars defines the OIFP.Curves noted with filled circles show annual mean MO dates.The black is the MO dates from the advanced horizontal range algorithm (AHRA) used in this study while red and orange are for early MO (EMO) and continuous MO (CMO) from the passive microwave algorithm (PMW) [6].Areas on the right axis give the extent of annual sea ice retreat within the region.

Figure 9 .
Figure 9. Mean annual evolution of the sea ice melt season for the Northern Hemisphere (March 1979-February 2017) for Case mask.Colors for the shaded bars define the mean SLIP, IIFP, SGIP periods.The full span of the shaded bars defines the OIFP.Curves noted with filled circles show annual mean MO dates.The black is the MO dates from the advanced horizontal range algorithm (AHRA) used in this study while red and orange are for early MO (EMO) and continuous MO (CMO) from the passive microwave algorithm (PMW) [6].Areas on the right axis give the extent of annual sea ice retreat within the region.

Table 1 .
Acronyms, definition and description of dates and periods.

Table 1 .
Acronyms, definition and description of dates and periods.

Table 2 .
The basic statistics of valid cell numbers for Case allv and Case core, linear decadal trends a (day/decade) and their margin of error (ME) b (day/decade).

Table 3 .
The North Pole hole attributes and the record periods used during the period of March 1979-February 2017.

Table 4 .
Basic statistical attributes of MO, linear decadal trends (day/decade) and their margin of error (ME) (day/decade) from different averaging procedures.

Table 5 .
Means Arctic dates and periods derived from different record mask cases.Only valid records that also have valid DOR are used for computing the means of individual dates or periods.b Case comm: Only records with all valid dates for the same ice year are used for computing the means.
a Case mask:

Table 6 .
Basic statistical attributes of Arctic dates and periods anomalies, their linear decadal trends (day/decade), and MEs (day/decade) of the trends.
* indicates statistical significance at the 99% confidence level.