Seasonal predictability of the minimum sea ice extent (SIE) in the Laptev Sea is investigated using winter coastal divergence as a predictor. From February to May, the new ice forming in wind-driven coastal polynyas grows to a thickness approximately equal to the climatological thickness loss due to summer thermodynamic processes. Estimating the area of sea ice that is preconditioned to melt enables seasonal predictability of the minimum SIE. Wintertime ice motion is quantified by seeding passive tracers along the coastlines and advecting them with the Lagrangian Ice Tracking System (LITS) forced with sea ice drifts from the Polar Pathfinder dataset for years 1992–2016. LITS-derived landfast ice estimates are comparable to those of the Russian Arctic and Antarctic Research Institute ice charts. Time series of the minimum SIE and coastal divergence show trends of −24.2% and +31.3% per decade, respectively. Statistically significant correlation (r = −0.63) between anomalies of coastal divergence and the following September SIE occurs for coastal divergence integrated from February to the beginning of May. Using the coastal divergence anomaly to predict the minimum SIE departure from the trend improves the explained variance by 21% compared to hindcasts based on persistence of the linear trend. Coastal divergence anomalies correlate with the winter mean Arctic Oscillation index (r = 0.69). LITS-derived areas of coastal divergence tend to underestimate the total area covered by thin ice in the CryoSat-2/SMOS (Soil Moisture and Ocean Salinity) thickness dataset, as suggested by a thermodynamic sea ice growth model.
Adaptation to a changing Arctic climate relies in part on our capacity to predict sea ice. Since the start of satellite monitoring of the Arctic in the late 1970s, observations have shown a decrease of the average September sea ice extent (SIE) at a rate of 13.3% per decade (Fetterer et al. 2017, updated daily). The retreat of the pack ice and the transition to a sea ice–free summer in the Arctic have important implications for marine ecosystems (Arrigo et al. 2008; Frey et al. 2015) and human activities, including navigation (Stephenson and Smith 2015; Pizzolato et al. 2016), tourism (Dawson et al. 2014; Johnston et al. 2017), and the exploitation of natural resources (Houseknecht et al. 2012a,b). This requires the development of skillful seasonal sea ice forecasting techniques.
In addition to its downward trend, the September sea ice also displays important interannual variability, which represents a considerable challenge for seasonal forecasting. In the Study of Environmental Arctic Change (SEARCH) Sea Ice Outlook 2008–13 (Stroeve et al. 2014), a collection of seasonal predictions of the pan-Arctic September SIE show that the forecasts are more accurate when the observed SIE is near the linear trend and lose accuracy when the observed SIE is anomalous. Furthermore, dynamical models have not provided higher forecast skill than simpler statistical methods. For instance, pan-Arctic September SIE hindcasts initialized in May or June from the four different dynamical models discussed in Blanchard-Wrigglesworth et al. (2015) display a root-mean-square error (RMSE) of 0.43 to 0.58 million km2, whereas the Williams et al. (2016) simpler statistical hindcast method based on winter sea ice export has a RMSE of 0.37 million km2. Understanding physical mechanisms that provide seasonal predictability of the interannual variability of the SIE is therefore crucial in the development of forecasting methods. Producing forecasts relevant to end users and stakeholders further requires the ability to generate regional forecasts. Some advances either from observation-based methods (Petty et al. 2017) or model-based methods (Germe et al. 2014; Day et al. 2014; Cheng et al. 2016; Bushuk et al. 2017b) show skill, but there remains a need for a better understanding of the mechanisms that provide regional predictability. The present paper focuses on winter preconditioning as a source of sea ice predictability in the Laptev Sea.
Winter dynamic preconditioning—in other words, the history of pack ice large-scale circulation during the previous winter—offers valuable information to estimate anomalies of the minimum SIE relative to the long-term trend. Nikolaeva and Sesterikov (1970) developed a forecasting technique for the Laptev Sea in which they track the timing of sea ice formation. They report that the new ice that starts forming in February grows to an average thickness of 1.15–1.45 m by the end of May, which is on the same order as the climatological thickness loss due to thermodynamic processes during summer. Building on this, they make a regional forecast of the summer ice edge from estimates of the amount of ice that survives (formed before February) or melts (formed after February) in their regional domain. This idea is supported by the model-based study of Chevallier and Salas-Mélia (2012), who report that the area covered by ice thicker than 0.9 to 1.5 m provides potential predictability of the August and September ice. Additional evidence of the winter preconditioning effect is obtained by observations of the sea ice export out of the Laptev Sea (Krumpen et al. 2013; Itkin and Krumpen 2017). Over the 1991–2011 period, Krumpen et al. (2013) report a negative correlation of r = −0.65 between the anomalies of winter ice flux and the following summer ice extent. Using a regional coupled ice–ocean model forced with realistic atmospheric reanalysis in the winter and a climatology in the summer (thus suppressing variability of ice loss due to summer processes), Itkin and Krumpen (2017) find a similar negative correlation. At the pan-Arctic scale, Williams et al. (2016) find from observations that advection of sea ice away from the Eurasian coastlines late in the winter and sea ice export through the Fram Strait are both well correlated with the minimum SIE (r = −0.58 and r = −0.72, respectively). Skillful predictions of the minimum SIE based on winter preconditioning are in contrast with the idea of a predictability barrier for SIE forecasts (Day et al. 2014; Bushuk et al. 2017b). Based on perfect-model forecast experiments across five different global climate models (GCMs) with initialization in January, May, and July, Day et al. (2014) show that sea ice forecasts initialized in May lose skill more rapidly than those initialized in January or July. Furthermore, they report that September ice extent in northern Arctic regions (including the Laptev Sea) is only predictable from July-initialized forecasts. Using the Geophysical Fluid Dynamics Laboratory (GFDL) seasonal prediction model for retrospective forecasts over the 1981–2015 period, Bushuk et al. (2017b) find that regional summer SIE forecasts initialized before May show little to no skill in the Laptev and East Siberian Seas (and no skill before June in the Beaufort Sea).
Seasonal predictability can be partly explained by reemergence of sea ice state anomalies (Blanchard-Wrigglesworth et al. 2011; Bushuk and Giannakis 2015). Winter preconditioning is associated with the reemergence of sea ice thickness anomalies from the winter to the following fall; for instance, thinner ice will lead to earlier melt, which then allows for more absorption of solar radiation early in the spring and further amplifies the melt. Bushuk et al. (2017a) find that sea ice thickness anomalies are enhanced during the summer melt season, driven by ice–albedo feedback. In turn, fall sea ice anomalies tend to recur the next spring, as lower sea ice concentration and extent in the summer and fall allow for more absorption of solar radiation in the upper ocean, which delays freeze onset and inhibits sea ice growth. Sea ice thickness anomalies and subsurface ocean temperature anomalies manifest memory on a seasonal scale.
Wintertime anomalies in large-scale atmospheric forcing are linked with the Arctic Oscillation (AO), which is defined as the first mode of variability of sea level pressure (SLP) in the Northern Hemisphere (Thompson and Wallace 1998). The AO loading pattern is zonally symmetric, with a high pressure center in the Arctic and two low pressure centers in the North Atlantic and North Pacific. A positive phase of the AO is characterized by a weaker than average polar high and more frequent offshore winds on the Eurasian coastline, associated with increased sea ice advection away from the Eurasian coasts, enhanced Transpolar Drift Stream (Rigor et al. 2002), and enhanced ice export through the Fram Strait (Williams et al. 2016). However, the relationship between Fram Strait export and AO index is not necessarily stationary in time (Jung and Hilmer 2001; Smedsrud et al. 2017). Note that the study by Jung and Hilmer (2001) uses the North Atlantic Oscillation (NAO) index—which is closely related to the AO—and finds that the correlation between NAO and Fram Strait export is zero before the 1980s and becomes 0.7 starting in the 1980s. In this study, we choose to focus on the large-scale atmospheric forcing rather than a regional vorticity-based index (e.g., Dmitrenko et al. 2008) because potential predictability of the minimum SIE using the AO index as a predictor has been demonstrated (e.g., Williams et al. 2016; Ogi et al. 2016; Park et al. 2018). Williams et al. (2016) report that the connection between the AO and the pan-Arctic minimum SIE is only present after the early 1990s and that the phase of the AO does not significantly affect variability of ice export in the Arctic Ocean peripheral seas in the earlier period. This is consistent with the decrease in the mean ice thickness in the western Arctic from the mid-1980s to the early 1990s associated with a persistent positive phase of the AO (Tucker et al. 2001), also analogous to the AO-driven abrupt reduction of area covered by multiyear ice over the same period (Rigor and Wallace 2004). The thinning of the pack ice increases responsiveness of ice drift to atmospheric forcing, resulting in positive trends in sea ice drift speed (Kwok et al. 2013; Tremblay et al. 2015). Since the phase of the AO is thought to have some (sub)seasonal predictability, it could potentially further increase the lead time for summer ice predictions (Stockdale et al. 2015; Kryjov and Min 2016; Smith et al. 2018). We ask how this relationship holds at the regional scale and investigate the link between the wintertime phase of the AO and coastal divergence in the Laptev Sea.
We define coastal divergence as advection of the pack ice away from the coastlines or the landfast ice edge (see section 3a). We develop a method based on the Lagrangian Ice Tracking System (DeRepentigny et al. 2016; Williams et al. 2016; Newton et al. 2017) to identify areas of coastal divergence based on the Polar Pathfinder observational sea ice drift dataset (Tschudi et al. 2016). We first test the method (as a proof of concept) in the Laptev Sea, where the mean pack ice drift is away from the coastlines. This work represents a bridge between previous studies that quantified sea ice outflow across the northern and eastern boundaries of the Laptev Sea (Krumpen et al. 2013) and studies that have used an Eulerian approach to investigate the location of Laptev Sea wind-driven polynyas and ice production within them (e.g., Willmes et al. 2010, 2011; Preußer et al. 2016). Preußer et al. (2016) report a significant correlation (r = 0.69) between ice production in polynyas and the area flux out of the Laptev Sea. However, as shown by Krumpen et al. (2013), most of the sea ice being exported out of the region comes from the western and central part of the Laptev Sea, with little contribution from ice formed in the polynyas. Our study provides a Lagrangian perspective on polynya activity in the Laptev Sea. That is, in addition to identifying the initial location of new ice formation in the coastal polynyas, we advect this new ice to identify its final position at the end of the winter. The coastal divergence area in the winter can thus be understood as a proxy for the distribution of thin ice at the start of the melt season. To test this assertion, we estimate the ice thickness within the area of coastal divergence using a one-dimensional thermodynamic model (section 3b). We also speculate that in other regions of the Arctic, the connection between coastal divergence and export from the peripheral shelf sea may be more ambiguous than in the well-bounded Laptev Sea. Therefore, an explicit quantification of coastal divergence as proposed in this work may be easier to generalize to other peripheral seas.
The paper is structured as follows. Section 2 presents the observational data used in this study: sea ice drift, concentration, thickness, landfast ice extent, 2-m dewpoint temperature, snow depth, and AO index. In section 3 we present the method developed to calculate areas of coastal divergence, based on Lagrangian integration of sea ice motion vectors. In addition, we describe a thermodynamic sea ice growth model used to calculate ice thickness in regions of coastal divergence. In section 4, we present a statistical seasonal forecast model for the minimum sea ice extent in the Laptev Sea based on coastal divergence and present the main results from the study. Conclusions of the study are summarized in section 5.
a. Sea ice drift
We use the National Snow and Ice Data Center (NSIDC) Polar Pathfinder (PPF) version 3 sea ice motion vectors (Tschudi et al. 2016) in the Lagrangian tracking model. The choice of the PPF dataset is motivated by its spatial and temporal completeness. We use weekly averaged velocity components of the sea ice drifts gridded on the Northern Hemisphere 25-km-resolution original Equal-Area Scalable Earth Grid (EASE-Grid). Using weekly averaged velocity is the shortest time averaging that filters random errors while retaining a sufficient temporal resolution (Tschudi et al. 2016). Merged sea ice drifts result from an optimal interpolation (Isaaks and Srivastava 1989) of drift vectors from buoy data, satellite remote sensing data (mainly from passive microwave imagery on which a maximum cross-correlation algorithm is applied), and free drift estimates derived from geostrophic wind reanalysis when no other data are available (see Table 1). In the Laptev Sea, fewer buoy drifts are available; the PPF ice motion vectors therefore rely more heavily on satellite-derived drift and free drift estimates. Although they use a different sea ice drift product, Rozman et al. (2011) perform a validation study in the Laptev Sea shelf region and show that drifts derived from passive microwave using a maximum cross-correlation algorithm are well correlated to in situ drift measurements made from moorings. At the pan-Arctic scale, PPF ice motion vectors obtained by an interpolation of all sources except buoy data yield a root-mean-square error of 3.4 cm s−1 when compared to buoy velocities (Tschudi et al. 2016). It is also reported in the technical documentation of the PPF dataset and in the detailed analysis of Sumata et al. (2014) that the most accurate ice motion vectors are obtained during winter; drifts retrieved during the summer period being more prone to errors due to surface melt and increased atmospheric water content affecting the identification of ice parcels from passive microwave remote sensing. Furthermore, the Scanning Multichannel Microwave Radiometer (SMM/R), which was in activity from 1978 to 1987, is more heavily biased toward lower drift speeds due to a lower temporal sampling (Tschudi et al. 2016). Taking into account the positive trends in sea ice drift speed observed in the Arctic (Kwok et al. 2013), the relative error on the merged ice motion vectors is less important in last 25 years of the record than in the earlier part. For all these reasons, we select the winter period from 1992 to 2016 for conducting the present study; 2016 is the last year for which PPF sea ice drifts are available over the full winter.
b. Sea ice concentration
We use the daily NOAA/NSIDC Climate Data Record passive microwave sea ice concentration version 3 (Meier et al. 2017; Peng et al. 2013). This dataset is based on satellite observations from the Special Sensor Microwave Imager (SSM/I) and Special Sensor Microwave Imager/Sounder (SSMIS) radiometers and provides sea ice concentration values on a 25-km-resolution polar stereographic grid. The accuracy of the sea ice concentration values during the winter months is of approximately 5% compared to coincident airborne data or other satellite remote sensing data with higher resolution (Meier et al. 2018). In the summer, concentration estimates tend to be less accurate. We produced weekly sea ice concentration averages from 1992 to 2016, gridded on the 25-km EASE-Grid. From this dataset, we construct a time series of the minimum SIE in the Laptev Sea. A grid cell is considered ice covered when its sea ice concentration is above 15%.
c. Landfast ice
We use the landfast ice extent from the Russian Arctic and Antarctic Research Institute (AARI) operational ice charts gridded on a 1.25-km-resolution EASE-Grid 2.0 (Selyuzhenok et al. 2015). The gridded dataset is available weekly from October 1999 to December 2013, with an exception between January and July 2002 due to a shortage of data. The gridded dataset was produced for a study of the landfast ice extent in the southeastern Laptev Sea and therefore its coverage is limited between 125° to 139°E (from the Lena Delta to the New Siberian Islands). The AARI ice charts are produced manually by experts who map sea ice conditions based on satellite, airborne, and ship-based observations, meteorological stations, and buoys. More details about the accuracy of AARI ice charts can be found in Selyuzhenok et al. (2015).
d. Sea ice thickness
We use the merged CryoSat-2 (CS2) and Soil Moisture and Ocean Salinity (SMOS) sea ice thickness dataset (Ricker et al. 2017b), hereafter referred to as CS2SMOS, available from October through mid-April, for years 2010 to present. SMOS sea ice thickness is retrieved from L-band brightness temperature measured by the Microwave Imaging Radiometer using Aperture Synthesis (MIRAS) instrument (Tian-Kunze et al. 2016). It is available daily at a 12.5-km resolution and is accurate for ice thinner than 1 m (Kaleschke et al. 2015). CS2 sea ice thickness is retrieved by sea ice freeboard measurements from the Synthetic Aperture Radar/Interferometric Radar Altimeter instrument (SIRAL) (Ricker et al. 2014). It is available weekly and monthly at a 25-km resolution and is accurate for ice thicker than 1 m. The advantage of the merged product is to make use of the different sensitivities of the two sea ice thickness retrieval methods and to cover the entire ice thickness range (Ricker et al. 2017a). Validation campaigns with airborne electromagnetic thickness measurements (AEM) showed that in the Beaufort Sea the CS2SMOS product underestimated the thickness of substantially deformed ice by up to 0.86 cm but was in good agreement with the AEM thicknesses over first-year ice regions, and in the Barents Sea, for ice with a low level of deformation, CS2SMOS had a bias of −25 cm compared to AEM measurements (Ricker et al. 2017a). We use the weekly temporal resolution dataset and regrid it on the 25-km EASE-Grid. As melt ponds start covering the pack ice in the spring, thickness retrieval from satellite measurement stops working. For that reason no data are available after mid-April.
e. Temperature at 2 m
The thermodynamic sea ice growth model (see section 3b) is forced with monthly-mean 2-m dewpoint temperature Td from the ECMWF interim reanalysis (ERA-Interim, herein ERA-I) dataset (Dee et al. 2011). Lindsay et al. (2014) report that from seven different reanalysis product, ERA-I provides the smallest bias compared to observations for near-surface air temperatures in the Arctic. We calculate monthly spatial averages over the Laptev Sea by taking the weighted average of all grid cells within a box delimited by 75°–85°N and 115°–130°E. We use dewpoint temperature as an approximation for the snow surface temperature following Raleigh et al. (2013). In the ERA-I dataset, over our region of interest, the 2-m dewpoint temperature is systematically 2°–3° lower than the 2-m air temperature, which is in accord with Huwald et al. (2005) measurements of the difference between skin temperature and 10-m air temperature during the Surface Heat Budget of the Arctic Ocean (SHEBA) mission.
f. Snow depth
The snow depth used in the thermodynamic sea ice growth model is specified from the Sever Soviet airborne expeditions (Shalina and Sandven 2018). Most of the Sever campaign measurements took place from March to May over the 1960s to 1980s. The entire Arctic Ocean was sampled during the multidecade campaign, with a special focus on the Russian marginal seas during the 1980s. Based on the Sever data, Shalina and Sandven (2018) report March–May mean snow depth over level ice in the Laptev Sea to be 9.8 cm, with a standard deviation of 3.1 cm. One of the key advantages of the Sever dataset is that the landings took place over both multiyear ice and first-year ice, and snow depth measurements are provided for various categories of ice or snow conditions (e.g., snow depth on level ice). This is a considerable advantage compared to the Soviet North Pole drifting stations that were mainly located in the central Arctic and established on multiyear ice. Data from the North Pole drifting station were used in the widely used Warren et al. (1999) snow depth climatology, thus making it not well suited for snow depth over first year ice and in the Russian marginal seas.
g. Arctic Oscillation index
We use winter mean [December–April (DJFMA)] AO index time series calculated from the monthly mean AO index distributed by the National Weather Service (National Weather Service/Climate Prediction Center 2005, updated daily).
a. Coastal divergence
We define coastal divergence as the area of sea ice exported from the coastlines, with units of km2. Note that the formal definition of divergence is (in s−1), where A (in km2) is the area of the Laptev Sea domain (Fig. 1), dA (in km2) is the area of coastal sea ice export, and dt (in s) is the duration of the tracking. In this paper, we refer to coastal divergence as dA only, the area of coastal sea ice export. Since A and dt are constants for any given time series that we present, they do not affect the correlations that we calculate.
1) Lagrangian Ice Tracking System
We use the Lagrangian Ice Tracking System (LITS; DeRepentigny et al. 2016; Williams et al. 2016; Newton et al. 2017) to track motion of sea ice over the course of a winter. Following Kimura et al. (2013), an array of tracers (one per EASE-Grid cell) is initialized at time t0 in the Laptev Sea and beyond its boundaries (see Fig. 1). This is done in order to capture the movement of sea ice entering the Laptev Sea from other peripheral seas (e.g., from the East Siberian Sea). The LITS is forced with weekly sea ice drifts from the PPF. Tracers are advected with a weekly time step from time t0, or Julian week X (WX), to time tf on the first week of May (W18); see Figs. 2a and 2b for an example of tracer advection using LITS. We test different initialization times t0, from the first week of December of the previous year (W–4) to the last week of April (W17), in order to identify the value of t0 that gives the highest correlation between late winter sea ice divergence and the minimum SIE. As tracers are advected in a Lagrangian manner, they can be located at any x and y position within a grid cell. DeRepentigny et al. (2016) report median and third quartile errors of 7% and 16%, respectively, by comparing LITS-derived trajectories to actual buoy trajectories using a weekly time step. Given a mean sea ice drift speed of about 5 cm s−1 in the Arctic Ocean, this represents median error of 2, 11, and 47 km for trackings of 1 week, 5 weeks, and 22 weeks respectively.
2) Distinction between landfast and drifting ice
Once the advection over the entire tracking period is complete (t0 to tf), a selection criterion is applied to distinguish between landfast ice tracers and northward drifting ice tracers (respectively red and blue points in Fig. 2). We choose a cutoff velocity of ucrit ≤ 0.5 cm s−1, applied as an average over all time steps, to discriminate between landfast ice and drifting ice. Landfast tracers are then returned to their initial position. In the southeastern Laptev Sea where a validation dataset is available, the resulting landfast ice is in good agreement with the AARI ice charts (see section 4b below).
The area of coastal divergence (Fig. 2c) is quantified by first constructing two polygons that represent the landfast and drifting ice edges at the end of the tracking period. These polygons are subtracted from the Laptev Sea regional domain and the resulting polygon (or set of polygons) is the coastal divergence area. To define the landfast ice polygon, we locate all landfast ice tracers within the Laptev Sea domain after they were returned to their initial position. To define the drifting ice polygon, we construct a binary 25-km EASE-Grid assigning a value of 1 to any grid cell that contains a drifting tracer; otherwise the grid cell is assigned a value of 0. We then take the edge of the largest continuous area as the drifting ice polygon. This process causes the rejection of a certain number of drifting tracers that become disconnected from the main pack when returned to the binary 25-km grid (green points in Fig. 2c). On average, 115 out of 1100 (~10%) drifting tracers within the Laptev Sea domain are discarded. The method used to define the drifting ice edge introduces an error of up to one grid cell width (25 km) multiplied by a typical length scale for the Laptev Sea (1500 km); that is, 0.03 million km2, in addition to a potential overestimation of the coastal divergence area up to 0.07 million km2 from rejected tracers (≈ 115 × 625 km2).
b. Sea ice thermodynamic model
Under divergent flow field conditions, the ice thickness gradient is expected to go from very thin ice near the landfast ice edge to ice of maximum thickness (hmax) on the northernmost edge of the area of coastal divergence. We calculate hmax using a steady-state one-dimensional two-layer (one layer of sea ice and one layer of snow) thermodynamic model with linear temperature profile in each layer (Semtner 1976). At the ice–ocean interface, the evolution of sea ice thickness (hi) is represented by
where the density of ice , the ocean temperature , the latent heat of fusion of ice , and the thermal conductivity of ice are assumed to be constant (see Table 2). Also, Ti is the temperature at the ice–snow interface, where continuity of the heat flux takes the form of
The model is integrated with a 1-day time step using a forward Euler explicit numerical scheme. The initial condition for the ice thickness is m. The snow surface temperature is set equal to the 2-m dewpoint temperature [Ts(t) = Td(t)], prescribed by the ERA-I 2-m monthly mean dewpoint temperature (section 2e). Snow depth increases linearly in time so that it starts at cm and reaches cm by the end of the integration [where 9.8 cm is the mean snow depth on level ice in the Laptev Sea reported by Shalina and Sandven (2018)].
The model is run for 70 days (10 weeks) starting in the first week of February. We choose February as the start date since ice formed after that date is best correlated with the following summer sea ice extent (see section 4a for details). The end date is set by the last available thickness field distributed in the CS2SMOS dataset (second week of April).
The simulated maximum ice thicknesses (hmax) for the winters of 2011–16 are reported in Table 3. We provide uncertainty on the simulated thickness with respect to variations of the snow depth. Shalina and Sandven (2018) report a standard deviation of 3.1 cm for the snow depth on level ice in the Laptev Sea. Removing or adding 3.1 cm of snow causes a ±5-cm difference in ice thickness. Another source of uncertainty comes from the snow deposition parameterization. If the snow is deposited as the maximum load that sea ice can sustain given hydrostatic balance, the prescribed snow depth of 9.8 cm is reached after ~1 week and the simulated sea ice thickness is on average 22 cm thinner. Finally, ice floe interactions, ocean heat fluxes, wind effect, and the contribution of synoptic events are not considered in our thermodynamic model but would contribute to the variability of the simulated ice thickness.
4. Results and discussion
a. Winter overview
The average December to May (W–4 to W18) area of coastal divergence is 0.45 million km2. As a reference, the total area of the Laptev Sea domain is 1.58 million km2; therefore, the total area of winter coastal divergence occupies approximately 29% of the domain. This compares to the winter sea ice export estimates by Alexandrov et al. (2000) of 0.48 million km2 over the 1979–95 period and by Krumpen et al. (2013) of 0.35 million km2 over the 1992–2011 period. Note that both studies consider the winter export from October to May.
Over the 1992–2016 period, we calculate the Pearson’s correlation coefficient of linearly detrended anomalies of coastal divergence and anomalies of the following September SIE minimum (green curve in Fig. 3). Given the relatively small size of our sample (25 elements), we estimate the uncertainty associated to the correlation coefficient by running a bootstrap analysis using 1000 resamples from the coastal divergence and minimum SIE time series. The standard deviation of the correlation coefficient from the bootstrap resamples is around . We study different initialization times and we observe a maximum correlation plateau when coastal divergence is integrated up until May (W18), from a start date between W3 (end of January) and W7 (middle of February); for example, starting on the first week of February (W5) gives a value of (p value = 7 × 10−4). Note that the length of integration is not constant for the different initialization times; for example, coastal divergence in late April (W16 to W18) is poorly correlated, but over such a short period of time it is hardly distinguishable from the noise. Results become meaningful for length of integration over 4 weeks.
This negative correlation suggests that years of stronger coastal divergence are associated with a lower SIE minimum. We argue that coastal divergence late in the winter leads to formation of new sea ice that does not grow to a sufficient thickness to survive the following summer melt, which brings additional confidence in the mechanism of winter preconditioning that was proposed in previous studies. Our study suggests that close to 40% (r2 ~ [−0.63]2) of the interannual variability of the minimum SIE anomalies can be explained by observed late winter sea ice dynamics. This is in agreement with Williams et al. (2016), who find a correlation coefficient of r = −0.58 between the pan-Arctic minimum SIE and a synthetic backtracked ice edge that serves as a proxy for late winter coastal divergence. Krumpen et al. (2013) find a correlation coefficient of r = −0.65 between the February to May anomalies of sea ice export in the Laptev Sea and the regional anomalies of summer sea ice. In a modeling study, Bushuk et al. (2017a) correlate the May sea ice thickness (SIT) to the September sea ice concentration on a gridpoint basis and find a maximum explained variance of roughly 40% in the regions with the largest September sea ice concentration variability (the Beaufort, Chukchi, East Siberian, and Laptev Seas, and north of the Kara and Barents Seas).
We also test the sensitivity of the correlation to different final times of integration of the coastal divergence (purple curve in Fig. 3). We perform a second advection experiment using a fixed 10-week period and move it through the winter, providing coastal divergence estimates between December to early February (W–4 to W6), and March to May (W9 to W18). It appears that early winter coastal divergence (December to early January) is poorly correlated with the minimum SIE, with the correlation starts being statistically significant after the end of January (W4) and improving as more of the late winter (end of March and April) is taken into account. The coastal divergence between April and the beginning of May itself is an essential, but not sufficient, contribution to the correlation. This confirms that the total sum of coastal divergence in the late winter, from late January/early February to early May, is the best seasonal predictor for the minimum SIE. In other words, the sum of all late-winter synoptic events (which drive coastal divergence) needs to be considered to provide the best estimation of the winter preconditioning effect on the pack ice.
Our results are in agreement with findings from Nikolaeva and Sesterikov (1970) and Krumpen et al. (2013), who note that, at the Laptev Sea boundaries, including advection that takes place before February reduces the correlation with the summer SIE. These results bring some additional elements to the discussion about the May barrier of predictability proposed in Day et al. (2014) and Bushuk et al. (2017b), who find little to no skill in regional and seasonal model forecasts initialized before May. Although we observe some relationship between coastal divergence from late January to April (W4 to W13), the strength of the correlation is moderate (r = −0.42). To maximize the correlation, coastal divergence must be considered until the first week of May, which echoes the 1 May forecasts presented by Day et al. (2014) and Bushuk et al. (2017b). Since seasonal forecasts represent an initial-value problem, the winter preconditioning effect not being maximized until May can contribute to the understanding of the barrier of predictability. However, there are other potential sources of explanation for this barrier of predictability within the models. Chevallier et al. (2013) raise an interesting point by assessing how closely the SIT structure in the CNRM-CM forecasts resembles the initialized state versus the uninitialized GCM state. They find that for May-initialized forecasts, two months into the forecast, the SIT structure approaches the uninitialized state of the GCM, which indicates that summer processes can progressively override the memory of the pack ice with respect to sea ice thickness. Despite that, several models have seen their forecast skill improve by initializing with observed or realistic sea ice thickness fields (Chevallier et al. 2013; Guemas et al. 2016; Blockley and Peterson 2018). However, these studies do not investigate September ice forecasts initialized before May, and further analyses of GCM diagnostics will be required in order to understand the mechanism responsible for the seasonal predictability barrier found in some models. Preliminary work suggests that the large-scale atmospheric circulation and its variability is decorrelated from coastal divergence in some GCMs due to biases in the location of the climatological Beaufort high and Transpolar Drift Stream (DeRepentigny et al. 2016).
b. Landfast ice area and late winter coastal divergence
We focus on coastal divergence integrated from the first week of February to the first week of May (W5 to W18). Each year in the time period analyzed, coastal divergence occurs east of Severnaya Zemlya, north of the New Siberian Islands, and in the region north of the Lena Delta (Fig. 4). The southernmost edge of coastal divergence areas represents the location of polynyas. The width of total coastal divergence areas has important interannual and spatial variability. It ranges from tens of kilometers to hundreds of kilometers for the W5 to W18 advection experiment (Fig. 4). Our results are in accord with Willmes et al. (2011) and Preußer et al. (2016) regarding the location of coastal polynyas in the Laptev Sea.
AARI landfast ice estimates for the first week of May in the southeastern part of the Laptev Sea (125°–139°E) are superimposed onto the maps of LITS-derived landfast ice area (Fig. 4) for years 2000–13. The landfast extent in that region of the Laptev Sea is typically stable from 1 February to 1 June according to Selyuzhenok et al. (2015). We calculate time series of the area of landfast ice from the AARI ice charts and from our LITS-based method. Both time series do not show a significant trend. The mean landfast ice extent from AARI and from our method are respectively AAARI = 0.132 million km2 and ALITS = 0.116 million km2. We calculate anomalies of the landfast ice area relative to the respective mean of each time series (Fig. 5). The standard deviations of the landfast extent are σAARI = 0.006 million km2 and σLITS = 0.011 million km2. Although the variability of the LITS-derived landfast ice is higher, it follows well the direction of the anomaly compared to the AARI dataset (Fig. 5). The correlation between anomalies of AAARI and ALITS is r = 0.66, significant to the 95% confidence level, which shows the feasibility of retrieving landfast ice extent from a Lagrangian approach.
A shortcoming of the method is that it does not distinguish between actual landfast ice and sea ice that experiences no net motion over the tracking period. In the Laptev Sea those two cases coincide most of the time due to the mean pack ice drift flowing away from the coastlines, which contributes to retrieving a realistic landfast ice edge to the first order, especially over the February to May period when the landfast region is stable. However, in some years, including 1996, 1997, and 1998, regions tagged as landfast extend from Severnaya Zemlya toward the central Laptev Sea (Fig. 4). Rather than being grounded or attached to the coastlines, parts of these regions experienced very little drift. In regions where there is also convergence occurring along the coastline (e.g., north of the Canadian Archipelago), or where the mean flow is parallel to the coastline (e.g., the Beaufort Sea), a higher-resolution drift field will be required. These details should be kept in mind for future developments of the method.
The mean area of coastal divergence integrated from February to May (W5 to W18) over the 1992–2016 period is 0.24 million km2, approximately 15% of the Laptev Sea domain. A positive trend of +31.3% per decade is present and statistically significant at the 95% confidence level (Fig. 6). This positive trend in coastal divergence is in accord with positive trends in ice production in the Laptev Sea polynyas reported by Preußer et al. (2016). The standard deviation of the linearly detrended coastal divergence is σCD = 0.08 million km2. Over the same period, the mean minimum SIE in the Laptev Sea domain is 0.80 million km2, with a negative trend of −24.2% per decade, statistically significant at the 95% level (Fig. 6). The standard deviation of the linearly detrended minimum SIE is σSIEm = 0.22 million km2. The LITS-derived amount of coastal divergence represents twice the amount of ice area export as reported by Krumpen et al. (2013): monthly averaged total sea ice area flux across the boundaries of the Laptev Sea for the 1992–2011 period is between 0.03 to 0.04 million km2 for February, March, and April (~0.12 million km2 over the three months). Errors in the calculation of coastal divergence area [see section 3a(2)] could account for up to 60% of the difference. We also hypothesize that part of this difference is due to within pack ice convergence in the Laptev Sea domain as sea ice drifts northward. This hypothesis could be tested using LITS with the output diagnostics of a fully coupled ice–ocean model.
The best linear fit between anomalies of SIE minimum and coastal divergence yields a slope of m = −1.7 (Fig. 7a), indicating that positive feedbacks (e.g., ice–albedo feedback) amplify the late winter coastal divergence by a factor of 1.7 during the melt season. This is in accord with results from Bushuk et al. (2017a), who demonstrate that the GFDL Climate Model displays a robust enhancement of sea ice thickness anomalies during the summer. At the pan-Artic scale, Williams et al. (2016) note a fivefold amplification between the area of winter preconditioning and the September SIE. The magnitude of the amplification factor can vary from one Arctic peripheral sea to another. The amplification factor will depend on sea ice thickness, as thicker ice will lose less area due to solar input in the mixed layer for a given open water fraction; cloudiness, which will affect the solar input in the mixed layer; and storminess (Day and Hodges 2018), which breaks the floes and increases the area available for lateral melt. Ice thickness in turn will depend on the mean ice circulation in the late winter. Our observations-based results indicate that this sea ice anomaly enhancement is a process that needs to be well represented in models in order to simulate correctly the interannual variability in the minimum SIE and the link with large-scale modes of variability in the atmospheric circulation. Narrowing down the uncertainty about the exact value of this amplification factor will require further study.
To investigate the link between regional coastal divergence and atmospheric forcing, we compare the February to May (W5 to W18) coastal divergence with the winter mean AO index (Fig. 7b). The correlation coefficient between the December to May AO index and the anomaly in coastal divergence is (Fig. 7b). This is in accord with the results of Williams et al. (2016), who show a significant correlation (r = 0.59) between anomalies of sea ice transport away from the Eurasian coastlines and the AO index. However, it contrasts with Krumpen et al. (2013), who only find a weak correlation between the winter sea ice area flux and the AO index. Our results indicate that a positive (negative) phase of the AO, associated with enhanced (weakened) Transpolar Drift Stream, is associated with positive (negative) anomalies of coastal divergence in the Laptev Sea.
c. Comparison with the CS2SMOS ice thickness
We investigate the connection between winter sea ice dynamics and thermodynamic growth of sea ice through inspection of ice thickness in the Laptev Sea. The sea ice model described in section 3b is used to estimate the maximum thickness of level ice (hmax) that grows from zero during the divergence period. This value is then used to contour a thickness isoline in the observational CS2SMOS dataset (Ricker et al. 2017b), which is then compared to the divergence area calculated using LITS, for the years covered by both datasets (2011–16).
We integrate coastal divergence from February to mid-April (W5 to W14) when the CS2SMOS SIT dataset is available. Estimated values of hmax range between 1.09 and 1.29 m, with a mean thickness of 1.15 m (Table 3). We expect climatological summer thermodynamic processes to melt sea ice up to that thickness range or thinner (Nikolaeva and Sesterikov 1970; Chevallier and Salas-Mélia 2012). Note that the CS2SMOS dataset provides the mean thickness within a grid cell, whereas the thermodynamic model gives an estimate of the level ice thickness. However, this should not represent a leading-order source of error in the comparison between the two datasets: along the Eurasian coastline the mean flow field is divergent, which causes sea ice to be mainly level ice.
As expected, regions of coastal divergence from LITS collocate with ice in the [0, hmax]-cm thickness range in the CS2SMOS dataset (Fig. 8). However, the LITS falls short of capturing the entire distribution of thin ice in CS2SMOS. The exact position of the northernmost edge of the coastal divergence area does not match closely the h = hmax isoline in CS2SMOS, except for 2011 when the two are in good agreement. The differences appear to be the largest in the central Laptev Sea where the northernmost extent of the coastal divergence area is south of the hmax isoline. The area of ice with thickness h ≤ hmax in the CS2SMOS dataset is therefore underestimated by the LITS (Table 3). The average areas over the 2011–16 period are 0.27 and 0.44 million km2 for LITS and CS2SMOS, respectively.
Error sources include uncertainties in the observational sea ice drifts (see section 2a), errors in the CS2SMOS merged product (see section 2d), and errors in the simulated sea ice thickness. For the errors in the sea ice drifts, preliminary analysis of the PPF dataset indicates a low bias in the satellite-derived drift speeds with respect to buoy data, which will contribute to smaller estimates of thin-ice area from LITS (since slow-moving ice would have more time to thicken in place). Reinterpolation of the PPF dataset correcting for the biases could contribute to minimizing this source of error and will be the object of future work. Uncertainties arising from the thermodynamic model and the CS2SMOS dataset are more difficult to address as they can potentially offset one another: a faster snow build-up leads to ice estimates (section 3b) that are 22 cm thinner, and the CS2SMOS dataset can underestimate ice with a low level of deformation by up to 25 cm (section 2d). If we consider only the uncertainty associated to the thermodynamic model, 22-cm-lower values of h = hmax actually bring LITS-derived and CS2SMOS areas of thin ice much closer, with mean values of 0.27 and 0.26 million km2, respectively. Last, another potential explanation for this discrepancy comes from the fact that our method is designed to detect direct coastal divergence leading to flaw polynya opening; that is, we are not attempting to resolve divergence within the pack ice [see section 3a(2) for details]. Yet the pack ice may experience divergence as it is advected offshore, resulting in opening of leads where new ice forms and thereby creating areas of thinner ice that are not captured by a Lagrangian approach using low-resolution drift vectors. In other Arctic peripheral seas (e.g., the Beaufort Sea), this divergence within the pack ice also contributes to preconditioning for the melt season (e.g., Hutchings and Perovich 2015).
Although this first analysis seems to indicate that the LITS sea ice drift–based approach to calculation of coastal divergence tends to underestimate areas of thin ice when compared to an observational sea ice thickness dataset, we acknowledge that there are very large uncertainties and that the record of 6 years is too short to draw robust conclusions. This first analysis sets the bases for future work that should address the ice thickness distribution at the end of the winter and its influence on the summer melt in a Lagrangian perspective. Using data from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) reanalysis (Zhang and Rothrock 2003) to estimate the average ice thickness lost to summer thermodynamic process in the current climate will constitute a first step.
d. September sea ice extent hindcast
We build a hindcast model for the 1992–2016 period based on late winter coastal divergence (Fig. 9). We evaluate the departure from the linear trend in SIE using the anomaly of coastal divergence integrated from February to May (W5 to W18). The departure from the trend is the coastal divergence anomaly multiplied by the slope of the best linear fit between anomalies of coastal divergence and anomalies of minimum SIE, using a leave-one-out approach (i.e., the best linear fit is computed by taking the target year out of the time series). The standard deviation of the hindcast error is σerr = 0.20 million km2 and the explained variance of the hindcast is r2 = 0.59. This represents 21% more explained variance compared to a hindcast based only on persistence of the linear trend (r2 = 0.38). The minimum sea ice extent in the Laptev Sea is primarily driven by the negative trend in sea ice associated with anthropogenic climate warming, yet including coastal divergence as a seasonal predictor in a statistical model significantly improves the hindcasts as it captures some of the interannual variability.
We investigate the predictability of the minimum sea ice extent (SIE) in the Laptev Sea on a seasonal time scale, using coastal divergence as a predictor. Coastal divergence occurs when sea ice is advected away from the coastlines by surface winds. The method builds on the idea that larger areas of coastal divergence during winter are associated with extended areas of sea ice too thin to survive the summer melt, which results in reduced SIE the subsequent September. Following Nikolaeva and Sesterikov (1970), we developed a method that uses observationally derived sea ice drift vectors to quantify the area of coastal divergence during winter for years 1992–2016. We use the Lagrangian Ice Tracking System (LITS) (DeRepentigny et al. 2016; Williams et al. 2016; Newton et al. 2017) forced with sea ice drifts from the Polar Pathfinder (PPF) version 3 dataset to advect sea ice tracers using a weekly time step until the first week of May, with start weeks ranging from the first week of December of the previous year to the last week of April. The advection of tracers in a Lagrangian manner naturally separates the drifting pack ice, the landfast ice, and the area of thinner ice in between. The LITS-derived landfast ice in the southeastern Laptev Sea is in good agreement with landfast ice estimates from the Russian Arctic and Antarctic Research Institute (AARI) ice charts (r = 0.66). The mean LITS-derived and AARI landfast ice area (and standard deviation) are respectively ALITS = 0.116 million km2 (σLITS = 0.011 million km2) and AAARI = 0.132 million km2 (σAARI = 0.006 million km2).
Large and statistically significant trends are present for time series of minimum SIE (−24.2% per decade) and time series of February to May coastal divergence (+31.3% per decade). The correlation between linearly detrended anomalies of coastal divergence and minimum SIE is maximal (r = −0.63) when integration of coastal divergence starts in between late January and mid-February (Julian week 3 to week 7) and ends on the first week of May (Julian week 18). The anomalies of February to May coastal divergence also correlate well with the winter mean phase of the Arctic Oscillation (r = 0.69). We produce hindcasts of the regional September minimum SIE for the 1992–2016 period by evaluating the departure from the trend using anomalies of February to May coastal divergence. The hindcast provides an explained variance of r2 = 0.59, which is a 21% improvement on the explained variance obtained from persistence of the linear trend. This study demonstrates that dynamical winter preconditioning of sea ice provides predictability of the minimum SIE on a seasonal time scale and on a regional spatial scale. This increases our confidence in the physical mechanism linking motion of sea ice during winter to sea ice extent in summer via preconditioning, as presented in Nikolaeva and Sesterikov (1970) and also investigated in Krumpen et al. (2013) and Williams et al. (2016) using different metrics.
Using a two-layer thermodynamic sea ice growth model, we estimate the level ice thickness in mid-April (hmax) of the new ice that started forming in the coastal polynyas during the first week of February. The ice grows on average to a thickness of 1.15 m for the 2011–16 period, which is thin enough to melt via climatological summer thermodynamic processes. Our results are in agreement with previous findings from Chevallier and Salas-Mélia (2012), who identified the area covered by ice thicker than 0.9–1.5 m as providing potential predictability of the September sea ice area up to 6 months in advance. The total area covered by ice thinner than hmax in the observational CryoSat-2/SMOS (CS2SMOS) ice thickness dataset appears to be underestimated by LITS. Preliminary analysis indicates that biases in the PPF drift speeds contribute to that error, in addition to errors arising from uncertainties in the thermodynamic model and in the CS2SMOS observations. Future work will include a reinterpolation of the PPF dataset correcting for the biases in drift speed.
Future research will include using the LITS to calculate the same diagnostic presented in this study with outputs from global climate models. We are particularly interested in investigating the concept of a predictability barrier.
This research was supported by the Forecasting Regional Arctic Sea Ice from a Month to Seasons (FRAMS) project, which is funded by the Marine Environmental Observation, Prediction and Response Network (MEOPAR). This project received additional support from the Canadian Sea Ice and Snow Evolution Network (CanSISE), which is funded by the Natural Science and Engineering Research Council (NSERC) Climate Change and Atmospheric Research Program. This research is a contribution to project ONR-N00014-11-1-0977 from the Office of Naval Research. Charles Brunette is grateful for academic and financial support from the Fonds de Recherche du Québec–Nature et Technologie (FRQNT), ArcTrain Canada, Québec-Océan, and the Institut Nordique du Québec. Robert Newton was supported by the National Science Foundation (NSF-PLR 15-04404). We thank Valeria Selyuzhenok for kindly sharing the gridded AARI landfast ice dataset. Original AARI ice charts were obtained through the World Meteorological Organisation Global Digital Sea Ice Data Bank Project. Geometrical operations on polygons were performed using the University of Manchester Generic Polygon Clipper library, designed by Alan Murta. The merging of CryoSat-2 and SMOS data was funded by the ESA project SMOS+Sea Ice. We thank Matthieu Chevallier and two anonymous reviewers for their thoughtful comments that significantly helped to improve the quality of the manuscript.