Here, we will use output from the Pan‐Arctic Ice Ocean Modeling and Assimilation System (PIOMAS, see Supporting information S1 for details) to document the spatiotemporal variability of the sea ice within the LIA. It has undergone extensive validation with a variety of pan‐arctic observations, including thickness data (Schweiger et al., 2011 ; Wang et al., 2016 ; Zhang & Rothrock, 2003 ). With respect to the LIA, there is evidence that PIOMAS is unable to represent the thickness gradients that exist as one moves offshore, and as a result, it underrepresents thickness in the region (Schweiger et al., 2011 ). However, the interannual variability in the PIOMAS sea ice thickness within the LIA is similar to that from CryoSat‐2 satellite observations (Ricker et al., 2014 ; Figure S1 ).

Apart from these general statements, little is known about the processes responsible for the spatiotemporal variability of sea ice within the LIA. Recently, large‐scale climate anomalies have been observed to impact sea ice within the LIA, including the early collapse or absence of the Lincoln Sea ice arch at the northern end of Nares Strait (Kwok et al., 2010 ; Moore & McNeil, 2018 ), the formation of a previously unrecognized polynya in the Wandel Sea during February 2018 (Ludwig et al., 2019 ; Moore et al., 2018b ), and the absence of the Beaufort High during the winter of 2017 (Moore et al., 2018a ). These results suggest that ice in the region is more dynamic than previously thought (Loewen & Michel, 2018 ).

For the purposes of this paper, we will refer to the LIA as the region of old and thick ice along the margin of the North American and Greenland Arctic Basins, excluding the channels within the Canadian Arctic Archipelago. In situ observations of ice thickness in this region are limited, with the most comprehensive data coming from declassified submarine under‐ice sonar profiles (Bourke & Garrett, 1987 ; Lindsay & Schweiger, 2015 ). This data indicated that during the period from the 1960s to the 1970s, mean ice thickness was 5–7 m (Bourke & Garrett, 1987 ). The presence of the particularly thick ice in the LIA was attributed to ice motion associated with the Beaufort Gyre and the Transpolar Drift (Bourke & Garrett, 1987 ).

The oldest and thickest sea ice in the arctic is found in an area that extends along the margin of the Arctic Ocean from the Wandel Sea (north of Greenland), southwestward to the Canadian Arctic Archipelago (Bourke & Garrett, 1987 ; Lindsay & Schweiger, 2015 ; Maslanik et al., 2011 ; Melling, 2002 ; Tilling et al., 2018 ). Climate models suggest that this region will be the last to lose its perennial ice cover (Laliberté et al., 2016 ; Sou & Flato, 2009 ), thus providing an important refuge for ice‐dependent species (Folger, 2018 ; WWF, 2018 ). As such, knowledge of the climate of this so‐called Last Ice Area (LIA) as well as the variability and changing nature of its sea ice are important inputs for decisions regarding the establishment of marine protected sanctuaries (Krajick, 2010 ; Loewen & Michel, 2018 ).

The continued loss of Arctic sea ice (Stroeve et al., 2012 ) is cause for concern both as a dramatic indicator of our changing climate (Vihma, 2014 ) as well as for its impacts on fragile regional ecosystems (Hinzman et al., 2005 ; Post et al., 2009 ). In addition to the well‐documented reduction in sea ice extent (Parkinson & Cavalieri, 2008 ), there are also trends toward a thinner (Schweiger et al., 2011 ), younger (Maslanik et al., 2011 ), and more mobile (Spreen et al., 2011 ) Arctic ice pack. These changes have led to stresses on the entire spectrum of ice‐dependent organisms from ice algae to polar bears (Lange et al., 2015 ; Post et al., 2013 ). Climate models, under various emission scenarios (Sou & Flato, 2009 ; Wang & Overland, 2012 ), suggest that these changes will continue, resulting in additional stresses on fragile ice‐dependent ecosystems (Barber et al., 2017 ; Post et al., 2013 ).

2 Results

Figure 1 shows the annual mean pan‐arctic fields from PIOMAS for the period 1979–2018. The thickest sea ice (thicknesses > 3 m) can be found in the LIA, with two distinct regions where the sea thickness exceeds 4 m, referred to here as the LIA East (LIA‐E) and the LIA West (LIA‐W; Figure 1a). Please refer to the Supporting information for details on the boundaries of these two regions. Sea ice age data (Maslanik et al., 2011, Figure 1b) show that these two regions also have the oldest sea ice in the arctic (age > 5 years).

Figure 1 Open in figure viewer PowerPoint 2011 Sea ice: (a) thickness (shading, m) and motion (vectors, cm/s); (b) age (years); (c) concentration (shading, %) and motion (vectors, cm/s); and (d) magnitude of sea ice motion (shading, cm/s). Results in (a), (c), and (d) are annual means from Pan‐Arctic Ice Ocean Modeling and Assimilation System, 1979–2018; while results in (b) are April 30 median values from Maslanik et al. (), 1984–2017. In (a), the 3 and 4 m thickness isocontour are shown in black and white, respectively. In (b), the white contours encloses areas where the ice age is 5 years or older. Please refer to the Supporting information for details on the coordinates of the 4 m isocontours.

In contrast to the thickness and age data, the annual mean sea ice concentration within the LIA is uniform and close to 100% (Figure 1c). The annual mean sea ice motion field (Figures 1a and 1c) indicates the presence of the anticyclonic Beaufort Gyre and the southeastward motion of the Transpolar Drift (Serreze & Barrett, 2010; Thorndike & Colony, 1982). Annual mean sea ice motion data derived from an independent blend of observations and model output (Tschudi et al., 2016) is consistent with the PIOMAS results (Figure S2). Sea ice speeds within the LIA are low, likely because of relatively thicker and more compact ice of higher mechanical strength (Figure 1d).

Figure 2 indicates that on annual time scales, variability in sea ice thickness between the LIA‐E and LIA‐W regions is not strongly correlated, providing support for the identification of these regions within the LIA. In particular, the decorrelation length scale (De Benedetti & Moore, 2017), here defined as the median distance over which variability in ice thickness is correlated at the 0.8 level for the locations within the LIA‐E and LIA‐W regions where the annual mean PIOMAS sea ice thickness is a maximum, indicated by the “+” markers in Figure 2 is ~200 and ~400 km, respectively. To put these distances in context, the distance between these two locations is ~1,300 km.

Figure 2 Open in figure viewer PowerPoint Spatial correlation maps of the annual mean sea ice thickness with time series of annual mean sea ice thickness in the: (a) LIA‐E and (b) LIA‐W regions. In each case, the 0.8 correlation contour is showing indicating the region where variability in the time series in the respective regions can explain approximately 64% of the variance in the sea ice thickness field. The decorrelation length scale is defined as the median radius from the center point, defined by the “+,” to this contour. LIA‐E = Last Ice Area East; LIA‐W = Last Ice Area West.

This result suggests that sea ice variability in the two regions of the LIA should be considered separately. This is done in Figure 3, which shows the time series of the annual mean sea ice thickness, concentration, and ice motion in the two regions. Also shown are the linear trends over the period 1979–2018 as well as a measure of the interannual variability, the detrended standard deviation. All trends are statistically significant at the 99th percentile confidence interval using a test that takes into account the reduced degrees of freedom arising from the temporal autocorrelation of the time series (Moore et al., 2015). Please refer to the Supporting information for details on this significance test. In both regions, sea ice has been thinning at a rate of ~0.4 m/decade, resulting in a thickness loss from the late 1970s of ~1.5 m, a value commensurate with the peak‐to‐peak interannual variability. Sea ice concentration and ice motion trends and variability differ for the two regions, with LIA‐W featuring values approximately twice as large as LIA‐E. The corresponding time series for the entire Arctic Ocean north of 70°N (Figure S3) indicate that both the LIA‐W and LIA‐E are losing ice mass at a rate twice that of the basin average, a result consistent with Bitz and Roe (2004). In contrast, the trend in sea ice concentration loss with the LIA is between 25% and 50% of the basin average. Furthermore, ice speeds within the LIA‐W are accelerating at a rate twice that of the basin average, while those in the LIA‐E are accelerating at ~50% of the basin average.

Figure 3 Open in figure viewer PowerPoint Time series of annual mean sea ice: (a) thickness (m), (b) concentration (%), and (c) motion (cm/s) in the LIA‐E region from Pan‐Arctic Ice Ocean Modeling and Assimilation System, 1979–2018. Time series of annual mean sea ice: (b) thickness (m), (d) concentration (%), and (f) motion (cm/s) in the LIA‐W region from Pan‐Arctic Ice Ocean Modeling and Assimilation System, 1979–2018. For each time series, the linear trend is also shown as a dashed line with one detrended standard deviation above and below the linear trend indicated by the dotted lines. LIA‐E = Last Ice Area East; LIA‐W = Last Ice Area West.

To assess temporal changes in sea ice characteristics, the annual cycles of sea ice thickness, concentration, and motion in the two regions within the LIA are shown in Figure 4 for the first, 1979–1998, and second, 1999–2018, halves of the period of interest as well as the difference between the two halves. Sea ice in the two regions tends to be thinner, lower in concentration, and more mobile in late summer and early fall (July–September) as compared to the late winter and early spring (March–May). Over the period of interest, there is a divergence in the timing of the extrema in sea ice thickness between the two regions (Figures 4a and 4b). In particular, the minimum in the LIA‐W is tending to occur earlier in the season. The reason for this behavior can be seen in the difference in the seasonal cycle of sea ice thickness between the two periods. For the LIA‐W, this difference is largest during August at 1.3 m. To put this value in context, the mean difference in thickness over all the months was 0.8 m, with a standard deviation of 0.2 m; thus the August difference exceeded the mean difference by over two standard deviations. This difference is more uniformly distributed throughout the year in the LIA‐E. In contrast to the ice thickness results, the reduction in sea ice concentration tends to be more focused during the summer period with minima occurring approximately 1 month earlier that that for sea ice thickness (Figures 4c and 4d) In addition, the timing of the maximum ice speed in both regions is occurring earlier in the season, with marked increase in sea ice mobility in the LIA‐W that is not seen in the LIA‐E (Figures 4e and 4f).

Figure 4 Open in figure viewer PowerPoint The annual cycle in sea ice: (a) thickness (m), (b) concentration (%), and (c) motion (cm/s) in the LIA‐E region from Pan‐Arctic Ice Ocean Modeling and Assimilation System. The annual cycle in sea ice: (b) thickness (m), (d) concentration (%), and (e) motion (cm/s) in the LIA‐W region from Pan‐Arctic Ice Ocean Modeling and Assimilation System. In all cases, the annual cycles over 1979–1998 and 1999–2018 are indicated by the blue solid and dashed curves, respectively.The difference between the two periods is indicated by the red curves. LIA‐E = Last Ice Area East; LIA‐W = Last Ice Area West.

As noted above, there is considerable interannual variability in the annual mean sea ice thickness in the LIA‐E and LIA‐W regions. To investigate the environmental conditions responsible for this variability, the detrended annual mean sea ice thickness time series in the two regions were used to generate spatial correlation maps with the detrended annual mean sea ice thickness and two components of the sea ice motion (Figure 5). For both regions, a bimodal signature in the sea ice thickness correlations is seen, with positive values in the LIA and negative values along the Siberian coast. For the LIA‐E region, the correlations show that ice thicknesses variability in this region is associated with transport across the Arctic, southeastward in years with high thicknesses that would tend to advect ice towards this region (Figure 5a). Similarly, high thickness in the LIA‐W is also associated with southeastward motion, which is associated with a cyclonic anomaly in this instance, which would also tend to advect ice towards this region (Figure 5b).

Figure 5 Open in figure viewer PowerPoint Correlation maps of the annual mean sea ice thickness (shading) and ice motion (vectors) with the detrended time series of annual ice thickness anomalies in the (a) LIA‐E and (b) LIA‐W regions. The respective regions of interest, LIA‐E and LIA‐W, are indicated by the “+”. LIA‐E = Last Ice Area East; LIA‐W = Last Ice Area West.

To further explore the meteorological controls on these spatially coherent structures, the correlation between the detrended ice thickness time series from the two regions and the principal components of the three leading empirical orthogonal functions (EOFs) of the annual mean sea‐level pressure field north of 60oN were computed (Overland & Wang, 2010) using the ERA‐I reanalysis (Dee et al., 2011). The leading EOF represents the arctic oscillation (AO) (Thompson & Wallace, 1998), while the second and third represent the arctic dipole (Wu et al., 2006) and the Barents oscillation (BO; Skeie, 2000), respectively. The AO has been shown to drive gyre‐like sea ice motion within the Arctic Ocean (Rigor et al., 2002),while ice motion associated with the Transpolar Drift and export through Fram Strait have been shown to be associated with both the arctic dipole and BO (Skeie, 2000; Wu et al., 2006). With respect to the LIA‐E, the correlation coefficients with the three principal components are 0.38, 0.20, and 0.24 and so the circulation anomalies associated with all three EOFs contribute to the interannual variability in sea ice thickness. In the LIA‐W, the three principal components have correlation coefficients of 0.57, 0.04, and 0.16, and so it is primarily the AO with a smaller contribution from the BO that contribute to the interannual variability of sea ice thickness in this region.