Backward and forward drift trajectories of sea ice in the northwestern Arctic Ocean in response to changing atmospheric circulation

SOA Key Laboratory for Polar Science, Polar Research Institute of China, Shanghai, China Chinese Antarctic Center of Surveying and Mapping, Wuhan University, Wuhan, China College of Earth Ocean and Atmospheric Sciences, Oregon State University, Corvallis, Oregon Department of integrated physical and ecological modeling and forecasting, National Oceanic and Atmospheric Administration Great Lakes Environmental Research Laboratory, Ann Arbor, Michigan


| INTRODUCTION
In the Arctic amplification of climate warming (Screen and Simmonds, 2010), the extent and thickness of Arctic sea ice have declined dramatically in recent decades (Parkinson and Comiso, 2013;Lindsay and Schweiger, 2015), which dominated the sea ice loss at the global scale although the Antarctic sea ice extent shows a contrasting trend to increase (Simmonds, 2015). The Arctic amplification is not uniform for the whole troposphere. The ratio of the Arctic nearsurface warming to that of the whole tropospheric column is highest in autumn (Screen and Simmonds, 2010). The time of onset of ice surface melt and refreezing, determined by satellite passive microwave measurements, have shown that the Arctic-wide melt season has lengthened at the rate of 5 days/decade during 1979-2013, primarily because of delayed autumn freeze-up (Stroeve et al., 2014). The lengthened melting season resulted in the increase in melt pond coverage over the ice, which exhibits a slight upward trend of 2.4% from 2000 to 2011 (Rösel and Kaleschke, 2012).
Sea ice moves in response to ocean currents, wind stress, and the Coriolis force, as well as the effects of sea surface slope and internal ice stress (Tremblay and Mysak, 1997). The Arctic winds generate an anticyclonic circulation system in the Canadian Basin, that is, the Beaufort Gyre (BG), with average clockwise ice motion in the western Arctic Ocean; and a cyclonic circulation with its centre in the Laptev Sea or Barents Sea (Proshutinsky and Johnson, 1997). The Transpolar Drift Stream (TDS) flows between these two circulation cells, transporting sea ice from the north of Siberia toward Fram Strait. The BG and TDS are not isolated and their boundary is ambiguous because of regulation by the atmospheric circulation patterns (Kwok et al., 2013). The Beaufort High (BH) forces the oceanic BG (Proshutinsky et al., 2009). Collapse of the BH with anomalous low sea level air pressure (SLP) could result in an anomalous reversal of the normally anticyclonic surface winds and ice motion in the western Arctic (Moore et al., 2018). The circulation pattern of the BG is closely related to the seasonality of synoptic processes (Asplin et al., 2009). The collapse of the BH was more ubiquitous during summer when the frequencies of Arctic cyclone increases in the central Arctic, in contrasting to the lower Arctic cyclone activity during winter (Simmonds et al., 2008). However, during recent winters, more cyclonic activity extended from Barents Sea into the central or the western Arctic Ocean (Rudeva and Simmonds, 2015), which resulted in the collapse of the BH and the reversal of the BG (Asplin et al., 2009;Moore et al., 2018;Petty, 2018). The cyclonic activities also would lead to more moist intrusions into the Arctic, which is projected to increase in a warming climate (Screen et al., 2018). During these intrusions the near-surface air temperature (SAT) in the region close to the North Pole has risen above freezing point even during winters (Graham et al., 2017). On the other hand, both the reversal of BG and the enhanced positive polarity of the Arctic Dipole Anomaly (DA) can accelerate the northward advection of sea ice from the BG region into the TDS region (Asplin et al., 2009;Wang et al., 2009). During the reversals of BG, the ice from the BG can be advected northeastward into the TDS (Lukovich and Barder, 2006); while the enhanced positive DA tends to transport the ice from the BG to the TDS through the western boundary of BG .
Increase in air temperature and melt pond coverage over the ice can effectively reduce the ice strength at the floe scale by increasing the brine and pore volumes within the ice cover (Timco and Johnston, 2002), while the drastic decline in ice thickness and concentration is hypothesized to weaken the ice strength at the basin scale (Hutter et al., 2018). The decreases in ice strength might result in enhanced Arctic sea ice mobility and increase its vulnerability to the storms. Storms coupled with increasing fetch because of the decreased ice concentration generate large waves which can propagate into the pack ice, thereby causing of the flexural swell and failure of sea ice (Asplin et al., 2012 and2014;Thomson and Rogers, 2014). This process fragments floes, which then more responsive to oceanic and atmospheric forcings and more susceptible to lateral melting. Numerical modelling of Arctic sea ice (Zhang et al., 2012) has shown that a decline in ice volume during the period 2007-2011 results in a 37% decrease in ice mechanical strength and a reduction of 31% in internal ice force, which in turn leads to increases in ice speed (13%) and deformation rates (17%) in comparison with the 1979-2006 mean, particularly in the western Arctic because there has been a rapid depletion of thicker, older ice. Zhang et al. (2012) stressed that, while model studies can shed light on the changes in the dynamic properties of Arctic sea ice, it is important to monitor the changes through satellite and in situ observations because the fundamental knowledge used for numerical modelling is largely based on a system dominated by thick, multi-year ice. Satellite data have confirmed that the spatially averaged drift speed within the Arctic Basin increased by 10.6 ± 0.9% decade −1 between 1992 and 2009 (Spreen et al., 2011). Sea ice motion is an important factor determining summer ice retreat in the Arctic Ocean due to redistribution of ice mass (Ogi et al., 2010;Kimura et al., 2013). Climate modelling results suggest that ice export through the Fram Strait can explain 28% of variance in long-term trends of sea ice thickness over the Arctic Ocean (Langehaug et al., 2013). Sea ice outflow through Fram Strait is regulated primarily by the east-west dipole pattern of SLP anomalies with centres of action located over the Barents Sea and Greenland (Tsukernik et al., 2010).
The most significant decline in summer Arctic sea ice extent has occurred in the sector stretching from Beaufort Sea to East Siberian Sea (Petty et al., 2016). However, sea ice area within this sector has substantial interannual variability attributable to the atmospheric circulation (Vihma et al., 2012;Petty et al., 2016). This sector of Arctic Ocean is the major investigation region of the Chinese National Arctic Research Expeditions in summers of 2008, 2010, 2014, and 2016. During these cruises, ice camps were established in the northwestern Arctic Ocean at 84.6 N and 144.3 W,87.3 N and 172.3 W,81.1 N and 157.4 W,and 82.8 N and 166.5 W,respectively (Figure 1). Considering the climatological mean ice drift, these sites are located in the boundary region between the BG and the TDS (Petty et al., 2016). To identify the source and destination of sea ice in this region can help to thoroughly understand the mass exchange of sea ice between the systems of BG and TDS, which may further influence sea ice loss in the entire Arctic Ocean. The drift trajectories of ice floes are crucial not only from a climatological perspective, but also provides context for determining pollutant migration over the Arctic Ocean. For example, the polymer composition of microplastics deposited in sea ice depends significantly on the drift paths of the sea ice (Peeken et al., 2018).
In this study, a satellite-derived product of daily ice motion vectors was used to identify the year to year changes in the backward and forward trajectories of ice drifting from the ice camps during 1979-2016 and their responses to changes in the atmospheric circulation pattern. The SAT derived from the reanalysis data and sea ice concentration derived from passive microwave remote sensing were used to explore the changes in climate and ice conditions along the trajectories, and their association with the variance of sea ice advection.

| DATA AND METHODS
At each ice camp, one sea Ice Mass balance Buoy (IMB; Richter-Menge et al., 2006) was deployed to track ice drift and to measure the seasonal evolution of ice mass balance. The IMB was equipped with a Navman Jupiter 32 global position system (GPS) receiver and it sampled positional data hourly with the horizontal accuracy better than ±10 m. The operational period of each IMB deployed in August 2008 and 2010 was approximately 11 months, while the operational period of the IMBs deployed in August 2014 and 2016 lasted longer than 1 year.
The Polar Pathfinder 25-km EASE-grid sea ice motion vectors dataset version 3 , available from the National Snow and Ice Data Center (NSIDC), provides Arctic sea ice drift vectors from October 1978 to February 2017. The sea ice drift vectors are obtained through calculating the correlation of groups of pixels between satellite image pairs and merging with both buoy drift vectors obtained from the International Arctic Buoy Program and the reanalyzed wind data from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR). The NSIDC ice drift vectors were used to reconstruct the backward and forward trajectories of the ice camp deployment sites for the period from August 1979 to August 2016. This ice motion product was selected to be used for the following reasons. (a) It allows identification of the long-term change of sea ice kinematics because it covers a more extensive time period than other products, for example, those from the Centre ERS d'Archivage et de Traitement (CERSAT; Girard-Ardhuin and Ezraty, 2012) and the Ocean and Sea Ice Satellite Application Facility (OSI-SAF; Lavergne and Eastwood, 2010  F I G U R E 1 Deployment sites of the CHINARE ice camps (triangles) and the trajectories of the buoys deployed at the camps. Also shown are the backward (to 31 August of the previous year) and forward (to 31 August of the following year or to the end of the buoy operation) trajectories estimated using the NSIDC product. Open and solid circles denote the monthly positions and terminal points of a trajectory, respectively. Red stars show the locations where the SLP was used to calculate the CAI. The solid yellow sector indicates the climatological domain of the BH defined by Petty et al. (2018) OSI-SAF products. (d) Its uncertainty is comparable with the OSI-SAF product and smaller than the CERSAT product (Sumata et al., 2014).
To reconstruct the forward trajectories, we projected the initial positions of the camps to the same Cartesian coordinate system as used by the NSIDC ice motion vectors, which is centred on the North Pole with the Y axis following the Greenwich meridian and the horizontal resolution of 25 km. The zonal and meridional coordinates of the initial camp site are X(0) and Y(0), with the ice drift vectors at time t along the ice trajectories are U(t) and V(t). The zonal and meridional coordinates of the reconstructed ice trajectories at the time t were calculated as: where the δt is the calculation time step (1 day). Finally, the reconstruct ice trajectory in the Cartesian coordinate was reprojected back to geographic coordinates. The procedure to reconstruct the backward trajectories was the same as that for the forward trajectories but using a negative time step. The terminal points of the reconstructed backward and forward trajectories were set on 31 August of the previous year and the following year relative to the time of deployment of the buoys, respectively. This cutoff date was defined to make the study period close to 1 year for both backward and forward trajectories. Due to data gaps in the NSIDC ice motion product, the ice trajectories in the ice season 1992/93 and 1993/94 cannot be reconstructed. Using the data from August 1979 to August 2016, we reconstructed 280 ice trajectories in total for four camp sites. Among these cases, there were seven for which the reconstructed forward trajectories did not reach the defined cutoff date because the trajectory reached the ice edge prior to the cutoff time. These trajectories were eliminated from further statistical analysis. Comparison of the actual buoy drift trajectories and estimated forward trajectories showed good agreement in terms of the general direction and destination ( Figure 1); hence, we have confidence in applying the NSIDC ice vectors to reconstruct the drift trajectories camps would have experienced in other years. Changes in atmosphere and sea ice conditions during 1979-2016 along the fixed backward and forward trajectories were obtained by interpolation of the ERA-Interim reanalysis (Dee et al., 2011) 2-m SAT and multichannel microwave radiometer -special sensor microwave imager sounders (SMMR-SSMIS) sea ice concentration retrieved using the National Aeronautics and Space Administration Team algorithms (Markus and Cavalieri, 2000) onto the trajectories. The fixed backward and forward trajectories were reconstructed in the year of buoy operation. Several methods have been used to define the ice melt season based on the ice surface condition (Stroeve et al., 2014), the ice concentration (Parkinson, 2014), or the SAT (Lei et al., 2018). In this study, to investigate the climate change, the onset dates of melt and freezing points in early summer and autumn were identified using the reanalysis SAT. Following Lei et al. (2018), the melt and freezing points were set as 0 C corresponding to the melt of snow or surface near-fresh ice and −1.8 C corresponding to the freeze-up of seawater, respectively. Due to its synoptic-scale variability, generally associated with the cyclonic activities (Simmonds et al., 2008;Simmonds and Keay, 2009), the SAT falls below −1.8 C in autumn and rises above 0 C in early summer several times. To account for this, the onset dates were defined as the dates when the running mean SAT in a given temprol window was over −1.8 or 0 C. Šmejkalová et al. (2016) used a 31-day sliding temporal window to characterize lake-ice phenology for the circumpolar Arctic regions. For time series used in this study, the window of 10 days is the shortest window which can avoid the running mean SAT passing across the threshold temperature repeatedly. Changing the length of the sliding window from 10 to 30 days leads to variance in the onset dates of melt or freezing season of less then 3 days, which implifies the onset date is not very sensitive to the window size. Hence, we used the modest length of 20 days for the sliding window to identify the onset dates of melt or freezing seasons. Long-term changes in average SAT in the freezing and melt seasons were then obtained. The freeze-up period in the Arctic is among the windiest and stormiest time of the year (Sotiropoulou et al., 2016). The strong winds would induce the wind-wave driven fractures of the new ice, creating pancakes. Consequently, the stormy conditions can delay the thickening of the new ice; however, this does not preclude the possibility of freeze up occurring during relatively quiescent periods between storms. Thus, the freezing onset date identified using SAT can indicate atmospheric surface cooling that preceds and occurs during freeze-up. However, the threshold of −1.8 C cannot indicate the onset of basal freezing for multiyear ice because of the thermal inertia exerted by snow and ice. Taking the 2010 camp as an example, the basal ice melt of the 1.49-m ice covered by 0.15-m of snow ceased by October 16, 2010 when the 20-day running mean SAT fell below −11 C (Lei et al., 2018). Thus, the isotherm date of −11 C, defined as the time when the 20-day running mean SAT fell below −11 C, was used as an estimate for the onset of basal ice freezing. Choosing an SAT threshold for freezeup is somewhat arbitrary because the onset of basal freezing also depends on the thicknesses of snow and ice and on the heat flux from the upper ocean. Changing the threshold SAT from −7 to −15 C leads to the variance of 24 days for the isotherm date. However, the long-term trends of the isotherm dates based on the thresholds from −7 to −15 C reveal the same pattern as found for the isotherm date using −11 C.
For example, along the forward trajectory from the 2008 camp site, the isotherm dates corresponding to the thresholds ranging from −7 to −15 C have the long-term trends within a narrow range from 2.8 to 3.2 days/decade during 1979−2016. This gives us further confidence to use a constant threshold of air temperature for all study sites. Moreover, it is convenient to identify the spatiotemporal change in atmosphere forcing for ice growth by using a constant threshold.
To quantify the change in ice conditions, this study focused on the September average ice concentration along the September leg of the backward trajectories. This was because the central Arctic Ocean was covered by compact sea ice with less interannual variability in ice concentration in other months. To explore the influence of sea ice advection on the atmospheric thermodynamic forcing and ice conditions along the trajectories, we also calculated the average SAT in the freezing season and the September ice concentration along the backward trajectories obtained for the years 1979-2016. We focused on the backward trajectories because the geographical locations of these trajectories can be depicted using the latitude to a great extent (Figure 2), and because the meridional variance of Arctic atmospheric forcing and ice conditions are generally larger than their zonal variance. The statistical relationships between average SAT in the freezing season or the average September ice concentration against the average latitude of the trajectories were then obtained. To quantify the regulating effect of the atmospheric circulation patterns on the sea ice motion, we calculated the seasonal Arctic Oscillation (AO) and DA indices, as well as the monthly Central Arctic Index (CAI) and BH anomaly. The SLP data of the NCEP/NCAR reanalysis (Kalnay et al., 1996) north of 70 N were used to derive the empirical orthogonal function modes, of which the AO and DA are the first and second modes ). The CAI, which was defined as the difference in SLP between 84 N, 90 W and 84 N, 90 E as shown in Figure 1, was calculated using the European Centre for Medium-Range Weather Forecasts ReAnalysis (ERA)-Interim reanalysis data according to Vihma et al. (2012). We used the CAI, but not the east-west SLP gradient between the Barents Sea and Greenland as in Tsukernik et al. (2010), because we wanted to characterize the overall magnitude of ice motion through the TDS, not just for the region across the Fram Strait. The BH anomaly was calculated as the spatial average SLP anomaly over the domain of 75 -85 N, 170 E-150 W relative to the 1979-2016 mean using the ERA-Interim reanalysis data. This domain, as shown in Figure 1, was defined according to Moore et al. (2018).
Long-term trends of atmospheric patterns and ice records during 1979-2016 were determined using linear regression. Correlation analysis was employed to test the statistical relationships between pairs of parameters. However, some variables, for example, the latitude of reconstructed trajectory, the SAT, and the ice concentration, are expected to have a long-term trend. It is hard to distinguish the physically dependence and the long-term covariation when performing the correlation analysis on the raw data. To address this issue we also recalculated the regressions using the detrended data if the long-term trend of the raw data was significant at the 95% level. The level of significance of the trends and correlations was evaluated using a two-tailed t test.

| Interannual changes in the reconstructed ice trajectories
In the years 1979-2016, ice floes at the 2008 camp site were advected mainly from the Makarov Basin along the TDS, with a few cases moving from the Canadian Basin along the western boundary of the BG (Figure 2a). The CAI can explain 39% (p < .001) of the interannual variance of latitude of the original positions, that is, the meridional displacement along the backward trajectories ( Figure 3a). According to the reconstructed backward trajectory, the floe was advected rapidly westward along the southern boundary of BG in September-November 2007. This was consistent with the fact that the westward transport of sea ice during summer-autumn in the southern BG region was relatively high in 2007 compared with the 1979-2016 mean because of anomalous anticyclonic winds in the BG region (Petty et al., 2016). In the southern BG region, the rapid westward transport of sea ice was sustained into winter and spring of the 2007/08 ice season. This led to a large amount of sea ice in this region being advected to the East Siberian Sea, which further resulted in rapid retreat of the sea ice in the Beaufort Sea during the 2008 melt season because of the lack of thick multi-year ice (Kimura et al., 2013). Close inspection of the estimated ice trajectories showed the terminal positions of the forward trajectories depended mainly on the direction of ice motion in winter (December-February), because the ice floes originated from the camp sites would generally drift to the intersection region between the northeast of BG system and the downstream TDS region by that time. In those ice seasons (22 cases) in which the ice originating from the 2008 camp site clearly drifted into the BG region, the average winter AO index had negative polarity (−0.13) and the average winter BH anomaly was positive (3.6 hPa). Both characteristics can be related to enhanced anticyclonic motion of the sea ice in the western Arctic (Vihma et al., 2012;Moore et al., 2018). In those ice seasons (11 cases within red circles in Figure 2a) in which the ice clearly drifted into the TDS region, the average winter AO index had positive polarity (0.68) and the average winter BH anomaly was negative (−6.5 hPa). Both of these characteristics can be related to the enhanced cyclonic motion of sea ice in the western Arctic. The winter AO index or the winter BH anomaly can explain 16% (p < .05) or 27% (p < .01), respectively, of the variance of longitude of the terminal positions of the forward trajectories, that is, the zonal ice advection (Figure 4a). An enhanced TDS could have promoted advection of the ice floe from the 2008 camp site into the TDS region. The increased values of the CAI (p < .01; Figure 3e) and DA (p < .05) might be related to the lower longitude of the terminal positions by regulating the zonal ice advection, meaning the ice floe was advected into the TDS regions.
The ice camp of 2010 was located near the central axis of the TDS (Figure 2b). Thus, the reconstructed backward and forward trajectories during 1979-2016 were influenced primarily by the intensity of the TDS. The floes originated mainly from the northern East Siberian Sea or the Makarov Basin. The CAI can explain 34% (p < .01) of the variance in latitude of the original positions ( Figure 3b). However, no significant relationship was identified between the DA index and the latitude of the original positions. When the TDS was strengthened or when the BH was abnormally negative, the ice floes from the camp site tended to drift into the TDS region toward the Fram Strait; otherwise, they would be trapped within a region close to the North Pole or drifted into the BG region. The CAI and DA values can explain 46% (p < .001; Figure 3f) and 14% (p < .05), respectively, of the variance in longitude of the terminal positions. The average CAI value was 6.0 hPa for those cases in which the ice clearly drifted into the TDS region (13 cases within the red rectangle in Figure 2b). This value is much larger than the 1979-2016 average (0.33 ± 0.25 hPa). The winter AO index or the winter BH anomaly can explain 15% (p < .05) or 24% (p < .01), respectively, of the variance of the terminal positions of the forward trajectories (Figure 4b).
In comparison with earlier years, the sites of the ice camps of 2014 and 2016 were deployed at relatively low latitudes, that is, 81.10 and 82.75 N, respectively (Figure 2c This indicates that ice floes from the sectors of the Chukchi and East Siberian seas have tended to be advected into the central Arctic Ocean and subsequently into Fram Strait in recent years. Furthermore, the frequently occurring positive polarity of DA during recent summers would be likely to lead to the anomalous advection of warm air from the south and the flow of relatively warm waters from the North Pacific into the Arctic Ocean, reduced cloudiness and enhanced insolation, in addition to increased sea ice export, which could result in reduced Arctic sea ice extent, especially in the Pacific sector of Arctic Ocean Overland et al., 2012;Lei et al., 2016). Therefore, during recent summers, the ice has retreated early in the Chukchi and East Siberian seas, which is conducive to the accessibility of Arctic Northeast Passage (Lei et al., 2015).

| Climate change along ice trajectories
From the SAT interpolated along the estimated trajectories during 1979-2016, it was found that the onset dates of melt and freeze did not show significant difference among the trajectories, apart from one exception for the forward trajectory from the 2010 site that occurred early on 6 June ( Figure 5 and Table 1). This was because the ice from this site was advected southward to the Fram Strait. The long-term average of the onset dates of the melt determined using the SAT for the backward (forward) trajectories was 10-12 June (6-12 June), while the long-term average of the onset dates of the freezing determined using the SAT for the forward trajectories was 26-28 August. They did not show any significant long-term trend for any ice camp site. This is different from the long-term trend of surface freeze-up identified from satellite passive microwave measurements, which was about 10 and 8 days/decade during 1979-2013 in the Chukchi and East Siberian seas, respectively (Stroeve et al., 2014). Thus, the significant delays in autumn surface freezeup were mainly attributed to the ice-albedo feedback, but not the changes in autumn SAT. The increase in melt pond coverage over the ice during the summer also would enhance the ice-albedo feedback and facilitates the delayed freeze-up of ice surface (Lei et al., 2016). For the open water, the strengthened cyclonic activities during the summer and early autumn (e.g., Simmonds and Keay, 2009) also would promote the delayed freeze-up through inducing ocean vertical mixing (Rainville and Woodgate, 2009). The insignificant spatiotemporal change in the onset dates of freezing and melt points determined using the SAT implies that the poleward gradient for surface atmospheric thermodynamic forcing is weak during summer, and that the loss of summer Arctic sea ice cannot generate significant positive feedback for warming of the surface atmosphere over the ice zone. The melting of sea ice during summer would consume heat from the atmosphere, resulting in the SAT remaining near 0 C, which restricts local warming for the Arctic ice zone in summer. This mechanism was termed as the effect of waterice bath by Overland (2009). However, we suggest that summer warming in the central Arctic would be promoted once the ice-free region extended further northward into the central Arctic, because the ice-free region has reasonably weak ability to regulate feedbacks related to air temperature increase compared with the ice zone.
Based on the onset dates of freeze and melt, we defined the melt and freezing seasons as persisting from 12 June to 28 August and from 1 September to 6 June, respectively. The mean SAT in the melt season ranged from 0.5 to 1.0 C and did not show significant spatial variance among the estimated trajectories. The spatial variance of the SAT in the freezing season was relatively larger compared with the melt season. Its value obtained from the forward trajectories ranged from −18.6 to −20.6 C, which was slightly lower than that obtained from the backward trajectories (−18.2 to −19.0 C), with the exception for the 2010 camp. This exception can be explained by the fact that the forward trajectories from the 2010 site were located at relatively low latitudes compared with the other sites. The isotherm dates of −11 C along the forward trajectories were earlier than along the backward trajectories for all cases by 5-8 days, which indicates that the SAT during early autumn in regions north of the Canadian Arctic Archipelago and Greenland was lower than in regions north of the Beaufort, Chukchi, and East Siberian seas. Furthermore, the long-term trends toward earlier isotherm dates of −11 C were more significant along the backward trajectories than along the forward trajectories ( Figure 6). From the linear regression, the −11 C isotherm dates were delayed by 15-23 days from 1979 to 2016 along the backward trajectories, that is, much more than along the forward trajectories (11-14 days). This implies the onset of basal growth for the sea ice with the same thickness at low latitudes might be delayed further compared with the high-latitude pack ice zone. This mechanism, introduced by the warming SAT, would be strengthened by the ice-albedo feedback (Perovich et al., 2007), the increased downward longwave radiation because of the and −11 C and rises above 0 C, respectively strengthened poleward transport of moisture and the enhanced local evaporation from the leads (Lee et al., 2017), the increased heat flux transferring from the thin ice or leads (Screen and Simmonds, 2010), and the enhanced waveinduced break-up of the floes (Asplin et al., 2012 and 2014; Thomson and Rogers, 2014). Thereby, following northward expansion of the Marginal Ice Zone (MIZ) due to loss of Arctic sea ice in summer (Strong and Rigor, 2013), delayed onset of basal ice freezing would be exacerbated for the pan Arctic. 1980 1985 1990 1995 2000 2005 2010 2015 18 Sep.
28 Oct. Date 198019851990199520002005201518 Sep. 8 Oct. 28 Oct. Date 198019851990199520002005201518 Sep. 8 Oct. 18 Oct. Date 1980198519901995200020052015 18 Sep. During the freezing season from 1 September to 6 June, SAT increased significantly during 1979-2016 along both backward and forward trajectories originating from all the deployment sites (Figure 7). The ERA-Interim reanalysis revealed that the trend of the SAT north of 70 N in 1989(Screen and Simmonds, 2010 was relative marked in autumn and winter (1.6 C per decade) compared to other seasons (0.5-0.9 C per decade). The near-surface warming amplification during autumn-winter is aided by strong low-level stability which limits vertical mixing through the troposphere (Overland, 2009) and increase in downward longwave radiation (Lee et al., 2017). Data obtained from land-based meteorological stations north of 70 N also showed that Arctic temperature amplification has been more remarkable in autumn and winter than in summer, with the rate of increase of SAT in autumn and winter being 1.9-2.5 times that in summer (Chylek et al., 2009). This highlights the importance of seasonal variability Arctic climate change. The long-term trends were more significant for backward trajectories than for forward trajectories. This indicates that air temperature warming during the freezing season in the Arctic Ocean was more remarkable at low latitudes than at high latitudes. In summer, low-latitude regions are more likely covered by low-concentration sea ice. The open water or thin sea ice in these regions tend to have a greater rate of ice growth during the freezing season. This can lead to increased heat flux from the ice-ocean system to the surface atmosphere and increased the evaporative moisture from the leads to the lower troposphere, both of which can produce enhanced local warming during the freezing season (Overland, 2009;Screen and Simmonds, 2010).
Along the estimated backward trajectories in various years, the average SAT in the freezing season was significantly dependent on the average latitude of the trajectories at the statistical level of 99.9% (Figure 8). Sea ice advected from lower latitudes experienced higher SATs during the freezing season. Because both the average SAT and the average latitude of the trajectories in the freezing season have a significant long-term trend at the 95% level, we recalculated the regressions using the detrended data. The detrended SAT in the freezing season can be explained by the detrended average latitude by 51, 48, 77, and 68% for 1980 1985 1990 1995 2000 2005  Although the correlation coefficients based on the detrended data were slightly smaller than those based on the raw data, they were still significant at the 99.9% level. This further indicates that the influence of the intensity of the TDS on atmospheric thermodynamic forcing for sea ice growth is physically meaningful.

| Change in ice concentration along ice trajectories
Along the September leg of the fixed backward trajectories, the linear regression (Figure 9) shows that September ice concentration decreased significantly from about 95% in 1979 to about 40% by 2016 for the 2008 and 2016 camp sites, and to about 20% for the 2014 camp site (p < .001). However, the corresponding decreased value was relatively small from 90 to 65% during 1979-2016 for the 2010 camp (p < .05) because of its relative high latitude of the backward trajectories compared to other camp sites. Sea ice originating from low latitudes has tended to experience lower-concentration ice conditions during early autumn in recent years, which, associated with the increased melt pond coverage over the ice (Rösel and Kaleschke, 2012), would have led to increased absorption of solar radiation by the upper ocean. We assumed constant values of the transmittance ratio of solar radiation for both ice cover (0.04) and open water (0.93) using typical observed values of multiyear ice (Nicolaus et al., 2012) and leads (Pegau and Paulson, 2001). Thus, a decrease in ice concentration from 95 to 20% could lead to an increase from 0.09 to 0.75 in the composite transmittance of solar radiation into the upper ocean. This would increase oceanic heat flux under the ice and result in delayed onset of basal ice freezing.
Along the backward trajectories reconstructed in various years, the September ice concentration along the September legs was found to be significantly dependent on the latitude of trajectory (Figure 10). This implies the influence of the TDS on sea ice thermodynamic processes. Specifically, an enhanced TDS would advect more ice from low latitudes, with most from the Chukchi and East Siberian seas (Figure 2), where ice conditions would be conducive to increased input of solar heat into the ice-ocean system. For the 2014 and 2016 camps, the September ice concentration obtained from the backward trajectories in the 2007/08 ice season deviated significantly (by more than two SDs of the regression error) from the linear regression of ice concentration against latitude. This was because the September legs of the backward trajectories in this year were located in the eastern BG region close to Banks Island (Figure 2c,d)  Also shown are linear fits (blue lines) and their correlation coefficients and confidence levels, as well as the average latitude of the September trajectories the ice concentration was generally larger than that at the same latitude in the Chukchi and East Siberian seas. The SD for the ice concentration obtained from the various years decreased strongly from south to north. This was because higher latitudes will in general have higher sea-ice concentrations, and hence the range of variability would be expected to be smaller. We also recalculated the regressions based on the detrended ice concentration and latitude in September. The squared correlation coefficients based on the detrended data were 0.78, 0.73, 0.80, and 0.73, for the camp sites of 2008, 2010, 2014, and 2016, respectively, all of which were very close to those obtained using the raw data and significant at the 99.9% level. This increased our confidence on the results revealed using the raw data. For example, the detrended September ice concentration in the 2007/08 ice season also deviated significantly from the regression by more than two SDs.

| Uncertainties of the reconstructed ice drift trajectories
To investigate the influence of the accuracy of remote sensing ice motion products on the estimation of sea ice trajectories, we compared the distances between the positions measured by the buoys deployed at the ice camps and drift estimated using the NSIDC and OSI-SAF remote sensing products . We note that our buoy data have not been assimilated in any ice motion product. The OSI-SAF ice motion vectors with a time span of 48 hr are estimated using a continuous maximum cross-correlation method on pairs of satellite images (Lavergne and Eastwood, 2010). Because the OSI-SAF data were available only for the freezing season from October to April, the comparison focuses on this period using the starting positions measured by the buoys on 1 October. Starting from the positions of the 2008 and 2010 camps on 1 October, the distance between the buoy position and the position estimated using the NSIDC data was found to be comparable with that between the buoy position and the position estimated using the OSI-SAF data until late January, both of which have values smaller than 70 km (Figure 11)  F I G U R E 1 1 Distance between camp position measured by the buoy and the position estimated using the NSIDC (red line) or OSI-SAF (blue line) remote sensing products in the operational year of the buoy, and the average (±SD) distance in 2005-2015 (black line and grey shading) between the positions estimated using the NSIDC and OSI-SAF remote sensing products during 2005-2015 because there were gaps in the OSI-SAF data when the floes were advected into the MIZ. By mid-February, the average distance increased to about 170 ± 150 km in relation to the position of this camp. Thus, the NSIDC data have weaker skill to reproduce ice motion in the downstream region of TDS compared with OSI-SAF data, especially in the region close to Fram Strait. However, the gaps in OSI-SAF data in the MIZ, attributable to low ice concentrations, also limit the applicability of these data to the reconstruction of ice trajectories. In contrast, for the 2014 and 2016 camps, the distance between the estimated positions and the buoy position, and that between the estimated positions obtained from the two remote sensed products remained <100 km. This implies that both remote sensed products can reproduce the BG motion better than drift in the downstream region of TDS, and that the differences between two remote sensed products are not spatially uniform. This is likely attributable to the lower ice speed and the larger ice concentration in the BG region in comparison with the region close to Fram Strait. Both factors would lead to decreased uncertainty of the remote sensed ice motion products (Hwang, 2013;Sumata et al., 2014). Nevertheless, we suggest that NSIDC data are adequate for use in determining the terminal positions of the trajectories for ice advected from the camp deployment sites, because the data uncertainty remained reasonable when the ice floe drifted into the intersection region between the northeast of BG system and the downstream of TDS. Consequentially, the general direction of ice motion (toward the BG region or the TDS region) can still be clearly identified.

| Representativeness of the original positions of the reconstructed ice drift trajectories
In this study, we set the original positions at four deployment sites of ice camps to reconstruct the backward and forward ice drift trajectories. It is convenient to validate the reconstructed ice trajectories using the buoy data. However, this also raises questions about the representativeness of the defined positions because of the randomness of the camp deployment. To assess the representativeness of the defined positions, we reconstructed the backward and forward ice trajectories from 1979 to 2016 for 16 sites in a region of about 900 km × 500 km surrounding the camp sites, which also in the boundary region between the BG and the TDS. Firstly, we compared the forward ice trajectories from 20 sites, including the camp sites and the 16 other test sites, in four ice seasons of 2002/03, 2003/04, 2007/08, and 2012/13, which were selected because of their abnormal (larger or smaller than one SD related to 1979-2016 average) BH and CAI. In ice seasons 2002/03 and 2007/08, the annual (September to follow August) CAI was abnormal high, with the values of 3.4 and 4.2 hPa, respectively. Such pattern of atmospheric circulation is conducive to the advection of ice floes from most test positions to the downstream TDS system, as shown in Figure 12a,c. Furthermore, the forward ice trajectories show a cyclonic (anticyclonic) trend in the ice season of 2002/03 (2007/08) because of the abnormal low (high) winter BH anomaly with a value of -10.8 hPa (9.0 hPa). Thus, there were no ice floes advected distinctly into the BG system in 2002/03, while three from the southeast corner of the study region advected into the BG system in 2007/08. In contrast, the annual CAI was abnormal low, with the values of -0.8 and 0.5 hPa in the ice seasons of 2003/04 and 2012/13, respectively. Thus, there were no ice floes from the study region advected into the TDS system during both ice seasons. Especially in 2012/2013, all the ice floes were advected southwestward because of the abnormal high (10.4 hPa) winter BH anomaly. The abnormal low CAI occurring in all seasons in 2012/13 only with a mild exception in spring 2013, which resulted in an extremely weak sea ice outflow from the Arctic Ocean to Fram Strait throughout this year and consequently a marked recovery for the Arctic sea ice extent in September 2013, with annual minimum value of 5.10 × 10 6 km 2 , which exceeded that in September 2012 by about 50% .
Secondly, we performed the correlation analysis between the annual CAI or the winter BH anomaly and the difference in latitude of the backward trajectories or in longitude of the forward trajectories. Here, we used the difference in latitude/longitude, but not the latitude/longitude of terminal points, because the 20 test sites had a relatively large range in the original latitude/longitude. From comparison among the 20 test sites, we found that, the squared correlation coefficients obtained from the camp sites were close to those obtained from other test sites, and all the correlations shown in Figure 13 were significantly at the 95% level. For all test sites, the explained degree (R 2 ) of the annual average CAI for the variance of the difference in latitude from the original position to the terminal position of the backward trajectories is significant at the 99.9% level, which is more significant than that for the variance of the difference in longitude of the forward trajectories. This indicates the influence of the TDS on the backward trajectories was stronger than on the forward trajectories. Furthermore, the explained degree of the variance of the difference in longitude of the forward trajectories by the annual average CAI was more significant than by the winter BH anomaly, especially for the sites in the north. This implies the forward trajectories were more affected by the strength of TDS than by the strength of winter BG. The above statistical conclusions are consistent for both the camp sites and the other test sites. Thus, we can infer that the regulation mechanism of atmospheric circulations on sea ice movement identified from the reconstructed backward or forward ice trajectories originated from four camps is also suitable for a wide boundary region between the BG and the TDS.

| CONCLUDING REMARKS
The NSIDC remote sensed ice motion product was used to reconstruct backward and forward drifting trajectories from four camps established in the summers of 2008, 2010, 2014, and 2016. These camps were located in a region where ice drift can oscillate between motion in the BG or the TDS. Based on comparisons with buoy measurements and the OSI-SAF remote sensed product, we found the NSIDC product has sufficient capability to estimate the original and terminal positions of the backward and forward trajectories, although it underestimates the southward speed of motion when ice is advected into Fram Strait.
From the reconstructed backward and forward ice trajectories during 1979-2016, the original positions of all the camp sites were found to show significant tendency of increasing advection from lower latitudes in recent years. The floes from all the camps were advected to lowerlongitude regions and became more involved in the downstream of TDS system in recent years. An enhanced TDS would advect more ice from the lower latitudes of the Laptev, East Siberian, and Chukchi seas into the central Arctic Ocean, promoting the retreat of sea ice during summer in these peripheral seas of the Arctic. This mechanism also would enhance the ice-based material transport from Arctic shelf regions to the Arctic basin region. Furthermore, an enhanced TDS would play a crucial role in sea ice loss for the entire Arctic Ocean via the export of sea ice from the central Arctic Ocean into the Greenland Sea or the Barents Sea, where the ice would melt even in winter.
The original positions of the backward ice trajectories and the terminal positions of the forward ice trajectories were significantly related to the atmospheric circulation F I G U R E 1 2 The reconstructed forward ice drift trajectories originated from 20 sites in a wide region centred at the camp sites during four ice seasons. The black and red dots denote the start and terminal positions of the ice trajectories, and the black triangles denote the deployment sites of four ice camps. The values of winter BH anomaly and annual CAI are shown in each panel pattern. Generally, the CAI has better skill than the DA index in explaining the original and terminal positions of the backward or forward trajectories. The squared correlation coefficient (R 2 ) between the seasonal DA and CAI for the period of 1979-2016 was in the range 0.46-0.86, all of which were significant at the 99.9% level. This indicates both indices are related to similar patterns in nature. However, the DA index cannot indicate the orientation of the Arctic SLP dipole. The CAI is calculated across the TDS region, thus it can depict an east-west dipole pattern of SLP perpendicular to the TDS and can be closer related to the ice motion along the TDS compared to the DA index. Compared with the AO index, the BH anomaly has better skill in explaining the terminal positions of the forward trajectories. This is because the AO index characterizes the pan-Arctic atmospheric circulation north of 70 N, whereas the BH anomaly characterizes the atmospheric circulation over the western Arctic Ocean, which is related closely to the curl of surface wind stress and sea ice motion in this region. Under the scenario of an abnormal high BH anomaly, ice floes from the camps would tend to become trapped within a region close to the North Pole or they would drift into the BG region, which would restrict the advection of ice toward the Fram Strait. Thereby, such a situation is not conducive to the loss of sea ice over the Arctic Ocean. However, an abnormal high BH anomaly is conducive to the retreat of summer sea ice for the lower-latitude region of western Arctic through inducing anomalous Ekman drift of sea ice from the Beaufort, Chukchi, and East Siberian seas toward the central Arctic (Ogi et al., 2008).
During the freezing season, Arctic climate warming and the poleward gradient of SAT were found more significant than in the melt season. The melting of Arctic sea ice can damp an atmospheric climate warming signal in summer. The delayed surface freeze-up identified by satellite passive microwave observation is mainly attributed to the ice-albedo feedback and the wind-wave driven ice fracture and ocean vertical mixing, but not the change in SAT during autumn because no significant delayed trend has been identified for the onset of freeze-up point based on the SAT. The onset of basal freezing of sea ice at low latitudes is likely to have greater delay compared with the high-latitude pack ice zone. This is because low-latitude oceans absorb more solar radiation and they experience stronger positive ice-albedo feedback and wind-derived waves. Science Foundation grants OPP1722729 and OPP1740768. J.W. was supported by NOAA CPO Office of Arctic Research. Three anonymous reviewers are greatly thanked for their constructive reviews.