An Intercomparison of Satellite Derived Arctic Sea Ice Motion Products

: Arctic sea ice motion information provides an important scientiﬁc basis for revealing the changing mechanism of Arctic sea ice and assessing the navigational safety of Arctic waterways. For now, many satellite derived Arctic sea ice motion products have been released but few studies have conducted comparisons of these products. In this study, eleven satellite sea ice motion products from the Ocean and Sea Ice Satellite Application Facility (OSI SAF), the National Snow and Ice Data Center (NSIDC), and the French Research Institute for the Exploitation of the Seas (Ifremer) were systematically evaluated and compared based on buoys from the International Arctic Buoy Program (IABP) and the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) over 2018–2020. The results show that the mean absolute errors (MAEs) of ice speed for these products are 1.15–2.26 km/d and the MAEs of ice motion angle are 14.93–23.19 ◦ . Among all products, Ifremer_AMSR2 achieves the best accuracy in terms of speed error, NSIDC_Pathﬁnder shows the lowest angle error and OSI-405-c_Merged performs best in sea-ice drift trajectory reconstruction. Moreover, season, region, data source, ice drift tracking algorithm, and time interval all inﬂuence the accuracy of these products: (1) The sea ice motion bias in the freezing season (1.04–1.96 km/d and 11.93–22.41 ◦ ) is smaller than that in the melting season (1.13–3.90 km/d and 14.41–27.41 ◦ ) for most of the products. (2) Most products perform worst in East Greenland, where ice movements are fast and complex. (3) The accuracies of the products derived from AMSR-2 remotely sensed data are better than those from other data sources. (4) The continuous maximum cross-correlation (CMCC) algorithm outperforms the maximum cross-correlation (MCC) algorithm in sea ice drift retrieval. (5) The MAEs of sea ice motion with longer time interval are relatively smaller. Overall, the results indicate that the eleven remote sensing Arctic sea ice drift products are of practical use for data assimilation and model validation if uncertainties are appropriately considered. Furthermore, this study provides some improvement directions for sea ice drift retrieval from satellite data. and and


Introduction
In recent decades, Arctic sea ice has been changing rapidly in the context of global warming and Arctic Amplification, with important implications for global climate [1][2][3][4][5][6]. Sea ice drift is an important kinetic parameter of sea ice change, which not only affects the heat and momentum transfer between the ocean and the atmosphere, but also has an impact on resource development and navigation in polar regions [7][8][9].
Since the early 1980s, satellite remote sensing technology has gradually become an important tool for sea ice drift monitoring in polar regions which is complementary to limited in situ observations. The basic principle of sea ice drift retrieval based on remotely sensed data is to track sea ice in sequential images over the same area. In recent years, the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI SAF), French Research  [11] A few studies have compared several sea ice motion products [26][27][28][29][30]. The authors of [26] evaluated sea ice drift products from OSI SAF, Ifremer, and the product generated from Advanced Synthetic Aperture Radar (ASAR) in the Laptev Sea using Acoustic Doppler Current Profilers (ADCPs) buoys from 2007 to 2008. It was found that the statistical correlation coefficients between different products and field data ranged from 0.56 to 0.86. Additionally, the correlation coefficient calculated from the product retrieved by AMSR-3 of 23 E data published by Ifremer was higher than that from the OSI SAF product using the same sensor data. The authors of [27] validated three low-resolution and three mediumresolution sea ice drift products from OSI SAF and Ifremer with high-precision Ice-Tethered Profiler (ITP) buoys during 2008-2010 throughout the whole Arctic region. It was found that the accuracy of the product retrieved by the CMCC algorithm was significantly better than that of the MCC algorithm in the sea ice region with low drift speed. Low-resolution sea ice drift products usually underestimated the sea ice drift speed, especially in the melting season. The accuracy of sea ice drift products retrieved by SAR data was the best, while there was no obvious difference in the accuracy between products retrieved by scatterometer and radiometer data. Furthermore, the error of sea ice drift products was closely related to their time interval. The authors of [28] compared 2002-2006 sea ice products from NSIDC, OSI SAF, Ifremer, and KIMURA using International Arctic Buoy Program (IABP) buoy data in the Arctic Ocean. The results showed that the NSIDC product and OSI SAF product were more accurate than the Ifremer product; this phenomenon was more significant in regions with lower sea ice concentrations and slower sea ice drift speeds. Furthermore, the uncertainty of the sea ice drift product would increase with increases of drift speed. The authors of [29] verified the NSIDC_Pathfinder and OSI-405-c_Merged products of 2014 and 2016 with ice drifters deployed by the Chinese National Arctic Research Expedition (CHINARE) in the western Arctic Ocean. It was found that NSIDC_Pathfinder tended to underestimate the sea ice velocity, while OSI-405-c_Merged tended to overestimate it. Compared with OSI-405-c_Merged, NSIDC_Pathfinder had relatively lower error and higher temporal and spatial resolution, and thus, was more suitable for estimating sea ice deformation. Taking the IABP buoy data from 2009 to 2017 as the verification set, the authors of [30] evaluated and compared some mainstream satellite sea ice products from OSI SAF, NSIDC, and Ifremer in the whole Arctic. The results showed that the product accuracy was related with the spatial resolution of its source data, ice tracking algorithm, and drift merging method.
It was found that most previous studies which compared products related to satellite derived sea ice drift only evaluated the bias of speed and ignored the error of ice motion angle and the ability of the product to reconstruct sea ice drift trajectory. Moreover, most studies only verified the whole Arctic or a specific region. The inconsistency of ice drifts in different regions was ignored. Furthermore, as some sea ice motion products have been updated recently (e.g., NSIDC published a new version of Pathfinder in 2019), the conclusions from previous comparative studies may not be applicable for new products.
The objective of this study is to evaluate and compare eleven satellite products from OSI SAF, NSIDC, and Ifremer with high-precision buoys of IABP and Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) over 2018-2020. Speed error was used for validation, as well as angle error and the ability to reconstruct sea ice drift trajectory. Furthermore, the influence of season, region, data source, ice drift tracking algorithm, and time interval on the accuracy of the products was analyzed.

Sea Ice Motion Products
Eleven sea ice motion products for 2018-2020 were validated in this study, i.e., four products from OSI SAF, two from NSIDC, and five from Ifremer. Table 2 shows basic information regarding the eleven products. Sections 2.1.1-2.1.3. provide some details for the products from different institutions. Among the eleven products, most use polar stereographic projection with different projection centers and central meridians. As projection would cause length and angle distortion, which would affect accuracy assessments on drift speed and direction, in this study, the coordinate systems of all products were unified to the WGS-84 geocentric coordinate system. The four OSI SAF sea ice drift products validated in this study are OSI-405-c_Merged, OSI-405-c_AMSR2, OSI-405-c_ASCAT, and OSI-407-a [10,15]. OSI-405-c_Merged, OSI-405-c_AMSR2, and OSI-405-c_ASCAT are the low resolution sea ice drift products of OSI SAF with 62.5 km spatial resolution and two-day temporal resolution. The start and end time of the OSI-405-c series are centered at 12:00 UTC. OSI-405-c_AMSR2 and OSI-405-c_ASCAT are single sensor products derived based on the CMCC method from AMSR-2 and ASCAT, respectively. OSI SAF also provides another single sensor product which is retrieved from the SSMIS sensor. However, the OSI-405-c_SSMIS product is not validated in this study because this data set has a gap in 2018. Due to surface melting and moist atmosphere, it is difficult to track sea ice drift during summer based on passive microwave radiometers and active microwave scatterometers. Therefore, the single sensor products of OSI SAF only provide winter sea ice drift vectors from October to April. OSI-405-c_Merged is the multisensor product derived by merging the sea ice drift vectors retrieved from AMSR-2, SSMIS, and Metop-B ASCAT. OSI-405-c_Merged provides the whole year data. Details of the merging algorithm can be found in [15]. OSI-407-a is the medium resolution product with 20 km spatial resolution and daily temporal resolution, derived based on the MCC method. Band 4 (Infra-Red) of the whole year and Band 2 (visible) of May to September from the Advanced Very High Resolution Radiometer (AVHRR) data are used for derivation of OSI-407-a. OSI-407-a updates two times one day with start and end time centered on 6:00 or 18:00 UTC before 11 October 2018. After that, this product updates four times per day with start and end time centered on 6:00, 12:00, 18:00, and 24:00 UTC. Due to the cloud effect on optical satellite images, the spatial coverage of OSI-407-a is limited.

NSIDC Sea Ice Motion Products
The two sea ice drift products from NSIDC validated in this study are "Polar Pathfinder Daily 25 km EASE-Grid Sea Ice Motion Vectors" product (referred to as NSIDC_Pathfinder hereafter) and "AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness Temperatures, Sea Ice Concentration, Motion & Snow Depth Polar Grids" product (referred to as NSIDC_AMSR2 hereafter) [16,18]. The latest version (version 4.1) of the NSIDC_Pathfinder product was used, which was updated on April 2019, including 25 km daily and weekly motion products. For 2018 to 2020, SSMIS data, IABP buoys, and NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) reanalysis forecasts were used to generate daily estimates. SSMIS data were used to retrieve the initial motion vectors based on the MCC algorithm. Then, NCEP/NCAR wind-derived and buoy-derived motions were used to refine the initial vectors using the optimal interpolation method. More details about the algorithm for pathfinder ice motion retrieval could be found in [18,24]. The weekly motion is an average of all daily motion daily vectors retrieved in that week. As buoys are sparse points without sequential records at specific locations to generate average vectors, the weekly motion product could not be validated with buoys, and thus, was not included in this study. NSIDC_AMSR2 is a daily sea ice motion product generated from AMSR-2 data based on the MCC algorithm. The grid of NSIDC_AMSR2 is 12.5 km × 12.5 km. However, sea ice motion vectors are recorded every five grids, therefore the spatial resolution of NSIDC_AMSR2 is 75 km.

Ifremer Sea Ice Motion Products
This study validated five Ifremer sea ice motion products, namely, Ifremer_AMSR2, Ifremer_Merged, Ifremer_SSMIS_H, Ifremer SSMIS_V, Ifremer_ASCAT [11][12][13][14]. Ifremer_ AMSR2 series have three 31.25 km datasets (i.e., Ifremer_AMSR2_H, Ifremer_AMSR2_V, and Ifremer_AMSR2) with two-, three-, and six-day lags. Ifremer_AMSR2_H, Ifremer_ AMSR2_V are derived based on the MCC method from horizontally and vertically polarized 89GHz channels, respectively. Ifremer_AMSR2 is generated by combining Ifre-mer_AMSR2_H and Ifremer_AMSR2_V motion vectors. In this study, Ifremer_AMSR2_H, Ifremer_AMSR2_V were not used. Ifremer_AMSR2 with two-day lag was used to compare with other products. Ifremer_AMSR2 with different time intervals were intercompared to evaluate the effect of temporal resolution on the accuracy of derived motions. Ifre-mer_SSMIS_H, Ifremer SSMIS_V, and Ifremer_ASCAT are 62.5 km products with threeand six-day lags, deduced based on the MCC method from horizontally polarized 91GHz channel, vertically polarized 91 GHz channel, and Metop-B ASCAT, respectively. In this study, only three-day datasets of these products were validated. It should be noted that Ifremer also released a sea ice motion product from Metop-A ASCAT for 2018-2020. Given that OSI SAF only uses Metop-B ASCAT, the Ifremer product from Metop-A ASCAT was not validated in this study to keep consistency in the use of sensors. Ifremer_Merged is a 62.5 km product with time intervals of 3, 6, and 30 days. Three-and six-day lagged Ifre-mer_Merged datasets were generated by merging Ifremer_SSMIS_H, Ifremer SSMIS_V, and Ifremer_ASCAT (Metop-B). Details of the merging algorithm can be found in [14]. In this study, Ifermer_Merged with three-day lag was compared with other products. Three-and six-day Ifermer_Merged were intercompared to discuss the impact of temporal resolution on product accuracy. Monthly drifts are estimated by summing five consecutive six-day lagged drifts of the merged product starting on the first day of each month. Therefore, the monthly product could not be validated by buoys, and thus, was not used in this study.

Buoy Data
IABP and MOSAiC buoy data were used in this study to evaluate the accuracy of the sea ice drift products. The IABP succeeded the Arctic Ocean Buoy Project, established in 1978 [31]. Its objective is to maintain a network of automatic data buoys in the entire Arctic Ocean to monitor sea level pressure, surface air temperature, and sea ice motion [32][33][34]. The IABP dataset provides daily latitude and longitude locations of buoys at 0:00 UTC and 12:00 UTC. MOSAiC was the first year-round expedition into the central Arctic, drifting with the sea ice from September 2019 to October 2020 and conducting comprehensive measurements on the atmosphere, ocean, and sea ice. During the expedition period, various types of buoys were deployed on the sea ice in the vicinity of the icebreaker, such as Ice-Tethered Profiler (ITP) and Surface Velocity Profiler (SVP). The temporal resolution of the buoy released by MOSAiC ranges from 10 min to 3 h. In this study, 647 buoys comprising 476 from IABP and 171 from MOSAiC with continuous time records (i.e., buoys with record interruption of more than one month were filtered out) from 2018 to 2020 were used as the verification data set.

Validation Methods
Two methods were used for validation: error estimation of sea ice motion vectors and Lagrangian trajectory reconstruction. When conducting the two validation methods, buoy records that derive ice speeds higher than 60 km/d were filtered out from validation datasets. This filtering step referred to the product evaluation conducted by NSIDC [18]. The reason for the high buoy-derived speed may be the errors of the positioning system [25,35,36]. Furthermore, Ifremer_SSMIS_H, Ifremer_SSMIS_V, and Ifremer_ASCAT are attached with quality flags. According to the user's manual [11,12], negative flag values indicate the cases that the drift should be dismissed. Therefore, drift vectors with negative flag values were discarded before validation. Moreover, as the start and end time of the OSI-405-c series are centered on 12:00 UTC, in this study, it was assumed that products from NSIDC and Ifremer were also centered on 12:00 UTC. Therefore, only the 12:00 UTC latitude and longitude location records of buoys were used for validation. OSI-405-a product has two or four datasets centered on different times of one day. The different datasets were merged into one dataset centered on 12:00 UTC for validation based on a weighted mean method that can be given by where L 0 is the location of one grid center, D(L 0 ) is the merged vector at grid L 0 , D i (L 0 ) is the vector at grid L 0 from the dataset i, n is the total number of vectors at grid L 0 from the two or four datasets of the OSI-407-a product, w i is the weight associated with the vector D(L 0 ), and T C is the start tracking time (UTC) of vector D i (L 0 ).

Error Estimation
In the error estimation method, sea ice motion vectors of the product were first bilinearly interpolated to the start positions of buoys. Then, the speeds and angles of buoy-Remote Sens. 2022, 14, 1261 7 of 23 derived vectors and corresponding interpolated vectors were quantitatively compared using three indicators, which can be given by where V MAE and θ MAE are the mean absolute error (MAE) of ice speed and ice motion angle, respectively. S P and S B are the displacements of the interpolated vector from the products and the buoy at the same start location, respectively. θ P and θ B are the angles of interpolated and buoy-derived vectors at the same start location, respectively. ∆θ is the angle between the interpolated vector and the buoy-derived vector. T represents the temporal resolution of the product. N is the number of interpolated vectors. In addition, the correlation between the drift speed of buoys and the drift speed of the interpolated products, which is expressed by the Pearson correlation coefficient r V . r V , ranged from −1 to 1. The closer r V is to 1, the stronger the positive correlation of ice drift speeds between products and buoys.
In this study, in order to evaluate the impact of season, region, data source, ice drift tracking algorithm, and time interval on the accuracy of the product, V MAE and θ MAE were estimated not only in the whole Arctic for 2018-2020, but were also calculated in different seasons and subregions. Furthermore, the V MAE and θ MAE of the products with different data sources, tracking algorithms, and time resolutions were compared separately. In seasonal variation analysis, errors of the freezing season (October to April) and the melting season (May to September) were compared. In order to make full use of data in three calendar years of 2018-2020, not only data of October 2018-April 2019 and October 2019-April 2020 were used to calculate errors in freezing season, but also data of January to April 2018 and October to December 2020 were employed. As the OSI-405-c_AMSR2 and OSI-405-c_ASCAT only provide data in the freezing season, they were only compared with other products in the freezing season. In spatial variation analysis, the Arctic was divided into eight subregions, i.e., Central Arctic (CA), Chukchi Beaufort Seas (CBS), Laptev/East Siberian Seas (LESS), Kara/Barents Seas (KBS), East Greenland (EG), Hudson/Baffin Bays (HBB), Canadian Arctic Archipelago (CAA) and Bering Sea (BS), referring to [37] (Figure 1). As buoy numbers are limited in HBB, CAA, and BS, sea ice drifts were separately evaluated in the other 5 subregions. In the evaluation of three factors (i.e., data source, tracking algorithm, and time interval) on the accuracy, the controlled experiment was used in which other two factors were identical among the compared products except for the factor being tested.

Trajectory Reconstruction
Lagrangian trajectories were reconstructed from products under the Universal Polar Stereographic North coordinate system, taking positions of buoys at the first days of trajectories as the start points. The vectors on start points or intermediate points of reconstructed trajectories were generated based on the bilinear interpolation method. The similarity between trajectories from buoys and products were quantitively evaluated with two indicators, which can be given by where D e and D c are Euclidean distance and cosine distance between endpoints of the product-derived and the buoy-derived trajectories, respectively. , and thus the range of cosine distance is [0, 2]. The closer the cosine distance is to 0, the smaller the angle between the two vectors, that is, the more similar trajectory estimated by the product is to the trajectory of the buoy.
were estimated not only in the whole Arctic for 2018-2020, but were also calculated in different seasons and subregions. Furthermore, the V MAE and θ MAE of the products with different data sources, tracking algorithms, and time resolutions were compared separately. In seasonal variation analysis, errors of the freezing season (October to April) and the melting season (May to September) were compared. In order to make full use of data in three calendar years of 2018-2020, not only data of October 2018-April 2019 and October 2019-April 2020 were used to calculate errors in freezing season, but also data of January to April 2018 and October to December 2020 were employed. As the OSI-405-c_AMSR2 and OSI-405-c_ASCAT only provide data in the freezing season, they were only compared with other products in the freezing season. In spatial variation analysis, the Arctic was divided into eight subregions, i.e., Central Arctic (CA), Chukchi Beaufort Seas (CBS), Laptev/East Siberian Seas (LESS), Kara/Barents Seas (KBS), East Greenland (EG), Hudson/Baffin Bays (HBB), Canadian Arctic Archipelago (CAA) and Bering Sea (BS), referring to [37] (Figure 1). As buoy numbers are limited in HBB, CAA, and BS, sea ice drifts were separately evaluated in the other 5 subregions. In the evaluation of three factors (i.e., data source, tracking algorithm, and time interval) on the accuracy, the controlled experiment was used in which other two factors were identical among the compared products except for the factor being tested.  Trajectory reconstruction requires products to have good performance in terms of spatial continuity. Most products showed poor spatial continuity in melting season [38]; therefore, in this study, trajectory reconstruction was only carried out in the frozen season (i.e., 1 October 2018-30 April 2019 and 1 October 2019-30 April 2020). Even in these periods, not all products had good spatial continuity. Therefore, the average daily spatial coverage areas (Formula (7)) of the products were calculated and compared. In this study, the products with large spatial coverage areas were used for trajectory reconstruction.
where S is the average daily product space coverage area, n i is the number of sea ice drift vectors in the product of day i, L is the spatial resolution of the product, and N is the total number of days for trajectory reconstruction. As the records of some buoys are discontinuous in the periods for trajectory reconstruction, only 23 buoys were used for validation. Figure 2 shows all buoy trajectories for trajectory reconstruction. These buoys are evenly distributed, which ensures trajectory reconstruction can be carried out in most areas of the Arctic.
where S is the average daily product space coverage area, n i is the number of sea ice drift vectors in the product of day i, L is the spatial resolution of the product, and N is the total number of days for trajectory reconstruction. As the records of some buoys are discontinuous in the periods for trajectory reconstruction, only 23 buoys were used for validation. Figure 2 shows all buoy trajectories for trajectory reconstruction. These buoys are evenly distributed, which ensures trajectory reconstruction can be carried out in most areas of the Arctic.

Overall Error Analysis
The nine annual sea ice drift products were compared in an overall error analysis. Table 3 shows the MAE results of sea ice drift speed and angle for the eleven annual products. Ifremer_AMSR2 with a two-day time interval and Ifremer_Merged with a threeday time interval were used for comparison. Table 3. Overall error estimation results of nine annual products.  Table 3 shows that the speed MAEs of these nine products were 1.15-2.26 km/d. Ifremer_AMSR2 had the lowest speed MAE, i.e., 1.15 km/d while NSIDC_AMSR2 had the highest speed MAE, i.e., 2.26 km/d. Furthermore, the accuracies of the Ifremer products were very close to those of OSI SAF, i.e., all in the range of 1.0-1.5 km/d. Overall, products from NSIDC achieved low accuracy, i.e., 1.85-2.26 km/d, which may have been due to the one-day time interval which they use. As the spatial resolution of the data source is commensurate with or larger than the daily displacement of sea ice, the image noise would strongly affect the accuracy of the retrieved daily sea ice motion. In the two NSIDC products, NSIDC_Pathfinder achieved the greatest accuracy, which also reflects its advantages of exploiting multiple data sources including SSMIS, reanalysis wind, and buoy data.
The angle MAEs of the nine annual products were 14.93-23.19 • . NSIDC_Pathifinder achieved the lowest angle MAE, i.e., 14.93 • , which may have been due to its satellitederived motion vectors being refined by NCEP/NCAR wind-derived and buoy-derived motions. OSI-405-c_Merged, which uses the CMCC algorithm, had the second-lowest angle MAE, i.e., 15.07 • . A possible reason for this is that CMCC is an improvement on MCC, which can retrieve displacement lengths smaller than the image pixel size and effectively suppress quantization noise. Of all products from Ifremer, Ifremer_AMSR2 had the best angle accuracy, i.e., 17.09 • ; the others had an angle MAEs of nearly 20 • , which may have been due to their different spatial resolutions of data sources and products. Compared with Ifremer_AMSR2, NSIDC_AMSR2, which is also derived using AMSR2, data had the worst accuracy in terms of angle. The reason for this may be that it has a shorter time interval and only uses the V-polarization band to extract ice motion. Figure 3 presents boxplots showing absolute errors of drift speed (a) and angle (b) for the nine annual products. According to the boxplots, the medians of speed error obtained from OSI SAF (i.e., OSI_M and OSI-a) and Ifremer (i.e., Ifrem_AM, Ifrem_M, Ifrem_S_H, Ifrem_S_V, and Ifrem_AS) were similar, and OSI-407-a (i.e., OSI-a) had the lowest median error of speed. However, from the upper quartile, maximum and average values, Ifremer_AMSR2 (i.e., Ifrem_AM) had the lowest speed error. The medians of angle error from all products were similar except for NSIDC_AMSR2 (i.e., NSIDC_A). OSI-405-c_Merged (i.e., OSI_M), NSIDC_Pathfinder (i.e., NSIDC_P) and Ifremer_AMSR2 (i.e., Ifrem_AM) achieved similar low error distributions in direction. However, the MAE of NSIDC_Pathfinder was lower than OSI-405-c_Merged and Ifremer_AMSR2, which indicated that it achieved the lowest angle error. The conclusions from Figure 3 agree with those drawn in Table 3.   Figure 4 shows density scatter plots between ice speeds derived from ice drift products. It can be seen that the correlation r values of OSI-407-a NSIDC_AMSR2 (i.e., 0.54) were relatively small. The correlation coefficients ucts were all above 0.9, and Ifremer_AMSR2 achieved the highest r (i.e., 0.9 ure shows, the two daily products from NSIDC achieved more high-speed (s than 30 km/d) ice motions than other products. The possible reason for thi high speeds are averaged with low speeds over long time intervals in othe represents the number of pairs matched between buoy and product, which the spatial coverage of the product to a certain extent. The two products values are NSIDC_Pathfinder (i.e., 87,259) and OSI-405-c_Merged (i.e., 7 demonstrated that the usage of multisource data can increase the spatial c trieved sea ice motions.  Figure 4 shows density scatter plots between ice speeds derived from buoys and sea ice drift products. It can be seen that the correlation r values of OSI-407-a (i.e., 0.83) and NSIDC_AMSR2 (i.e., 0.54) were relatively small. The correlation coefficients of other products were all above 0.9, and Ifremer_AMSR2 achieved the highest r (i.e., 0.95). As the figure shows, the two daily products from NSIDC achieved more high-speed (speeds greater than 30 km/d) ice motions than other products. The possible reason for this may be that high speeds are averaged with low speeds over long time intervals in other products. N represents the number of pairs matched between buoy and product, which can represent the spatial coverage of the product to a certain extent. The two products with larger N values are NSIDC_Pathfinder (i.e., 87,259) and OSI-405-c_Merged (i.e., 78,873), which demonstrated that the usage of multisource data can increase the spatial coverage of retrieved sea ice motions. Figure 5 shows scatter plots between buoy-derived ice speed and absolute error of sea ice speed from each product. In this figure, only the buoys with speeds below 25 km/d were used and divided into five equal intervals (as shown in Figure 4) that best matched buoys move below 25 km/d. From Figure 5, it may be seen that the speed error of all products increased with the increase of ice speed. Furthermore, there was a strong positive correlation between the speed error and speed, i.e., all r (correlation coefficient) values were higher than 0.86. With the increase of ice speed, Ifremer-AMSR2 showed the smallest error increase, which reflects its superior accuracy.  Figure 5 shows scatter plots between buoy-derived ice speed and absolute error of sea ice speed from each product. In this figure, only the buoys with speeds below 25 km/d were used and divided into five equal intervals (as shown in Figure 4) that best matched buoys move below 25 km/d. From Figure 5, it may be seen that the speed error of all products increased with the increase of ice speed. Furthermore, there was a strong positive correlation between the speed error and speed, i.e., all r (correlation coefficient) values were higher than 0.86. With the increase of ice speed, Ifremer-AMSR2 showed the smallest error increase, which reflects its superior accuracy.  Figure 6 shows the density scatter plots between buoy-based ice speed and angle absolute error of each product. It was found that among all products, the large angle errors were mostly gathered in the very low speed regime. Referring to [27], the angle errors after removing slow buoys (i.e., speed lower than 3 km/d) were recalculated; the results  Figure 6 shows the density scatter plots between buoy-based ice speed and angle absolute error of each product. It was found that among all products, the large angle errors were mostly gathered in the very low speed regime. Referring to [27], the angle errors after removing slow buoys (i.e., speed lower than 3 km/d) were recalculated; the results are shown in Table 2. It can be seen that the angle errors of all products were reduced after deleting the buoy records with low speed. Moreover, the angle MAEs of NSIDC_AMSR2 and all products from Ifremer were greatly reduced (about 10 • ), but the decrease of the angle MAEs of OSI-405-c_Merged, OSI-407-a, and NSIDC_Pathfinder was small (lower than 5 • ), which may be due to the following reasons: (1) OSI-405-c_Merged adopts CMCC method which can extract sea ice drift at subpixel level, and thus performs better at areas with low ice speed; (2) The data source used by OSI-407-a is AVHRR with a resolution of 1.1 km, i.e., much higher than that of data sources used by other products, so the angle of low speed drift can be extracted more accurately; and (3) NSIDC_Pathfinder corrects the angle error with buoys and reanalysis data.  Figure 7 show the average daily spatial coverage area results of eleven products. It may be seen that the NSIDC_Pathfinder achieved the largest average daily spatial coverage values in both freezing and melting seasons, which demonstrates its superior spatial continuity. OSI-407-a showed the worst spatial continuity as it uses optical data which are highly affected by clouds. Most products showed poor spatial continuity in the melting season. For trajectory construction in freezing season, the six products with highest daily spatial coverage values were used in this study, i.e., OSI-405-c_Merged, OSI-405-c_AMSR2, OSI-405-c_ASCAT, NSIDC_Pathfinder, Ifremer_AMSR2, Ifremer_Merged. It was found that most trajectories could not be reconstructed successfully by the other five products.  Figure 7 show the average daily spatial coverage area results of eleven products. It may be seen that the NSIDC_Pathfinder achieved the largest average daily spatial coverage values in both freezing and melting seasons, which demonstrates its superior spatial continuity. OSI-407-a showed the worst spatial continuity as it uses optical data which are highly affected by clouds. Most products showed poor spatial continuity in the melting season. For trajectory construction in freezing season, the six products with highest daily spatial coverage values were used in this study, i.e., OSI-405-c_Merged, OSI-405-c_AMSR2, OSI-405-c_ASCAT, NSIDC_Pathfinder, Ifremer_AMSR2, Ifremer_Merged. It was found that most trajectories could not be reconstructed successfully by the other five products. Table 4. Average daily spatial coverage of each product (×10 6 km 2 ).

Product
The   Figure 8 shows an example result of trajectory reconstruc and Figure 9 shows the time series of daily Euclidean distance tance (blue line) for reconstructed trajectories of the six product ple, all products showed good performance in trajectory recons terms of Euclidean distance, the accuracy of the reconstructe lower than that of other products. The accuracies of reconstructe in terms of cosine distance. However, the Euclidean distances NSIDC_Pathfinder were slightly larger than those of other prod  Figure 8 shows an example result of trajectory reconstruction for the six products, and Figure 9 shows the time series of daily Euclidean distance (red line) and cosine distance (blue line) for reconstructed trajectories of the six products in Figure 8. In this example, all products showed good performance in trajectory reconstruction from Figure 8. In terms of Euclidean distance, the accuracy of the reconstructed trajectory was slightly lower than that of other products. The accuracies of reconstructed trajectories were similar in terms of cosine distance. However, the Euclidean distances of OSI-405-c-ASCAT and NSIDC_Pathfinder were slightly larger than those of other products. and Figure 9 shows the time series of daily Euclidean distance (red line) and cosine distance (blue line) for reconstructed trajectories of the six products in Figure 8. In this example, all products showed good performance in trajectory reconstruction from Figure 8. In terms of Euclidean distance, the accuracy of the reconstructed trajectory was slightly lower than that of other products. The accuracies of reconstructed trajectories were similar in terms of cosine distance. However, the Euclidean distances of OSI-405-c-ASCAT and NSIDC_Pathfinder were slightly larger than those of other products.   Table 5 shows the mean Euclidean distances and mean cosine distances of the reconstructed trajectories from the six products in two freezing seasons. It is found that OSI-405-c_Merged obtained the lowest Euclidean distances (i.e., 10.6 km and 37.3 km) in both freezing seasons. Furthermore, OSI-405-c_Merged obtained the lowest cosine distance (i.e., 1.5 × 10 −4 ) in the 2018-2019 freezing season and comparable cosine distance (i.e., 9.7 × 10 −4 ) with the lowest one (i.e., 6.7 × 10 −4 ) in the 2019-2020 freezing season. Therefore, OSI-405-c_Merged showed the best ability in trajectory reconstruction among these six prod-  Table 5 shows the mean Euclidean distances and mean cosine distances of the reconstructed trajectories from the six products in two freezing seasons. It is found that OSI-405-c_Merged obtained the lowest Euclidean distances (i.e., 10.6 km and 37.3 km) in both freezing seasons. Furthermore, OSI-405-c_Merged obtained the lowest cosine distance (i.e., 1.5 × 10 −4 ) in the 2018-2019 freezing season and comparable cosine distance (i.e., 9.7 × 10 −4 ) with the lowest one (i.e., 6.7 × 10 −4 ) in the 2019-2020 freezing season. Therefore, OSI-405-c_Merged showed the best ability in trajectory reconstruction among these six products. The reason for this may be that the accuracy of the reconstructed trajectory showed the combined error in speed and angle. From Table 3 and Figure 3 showing error distributions of the products, ifremer_merged achieved the lowest speed error, but its angle error was larger than that of OSI-405-c_merged and NSIDC_Pathfinder. For angle error, NSIDC_Pathfinder performed the best, but its speed error was large. OSI-405-c_merged had relatively low values in both speed and angle error, and thus, achieved good performance in trajectory reconstruction.  Figures 10 and 11 show the spatial distribution of the Euclidean distance and cosine distance between the reconstructed trajectory of each product and the buoy trajectory in two freezing seasons, respectively. From Figures 10 and 11, it may be seen that almost all reconstructed trajectories of OSI-405-c_Merged and OSI-405-c_AMSR2 obtained low Euclidean distances and cosine distances from buoy trajectories, which demonstrated their good ability in trajectory reconstruction. This conclusion agrees with those drawn in Table 5. Furthermore, from Figure 10, it may be seen that the Euclidean distances of almost all products between the East Siberian Sea and the North Pole were larger than those near the Beaufort Sea and in the north of Ellesmere Island. However, from Figure 11, it may be seen that the cosine distances located between the East Siberian Sea and the North Pole were smaller than those near the Beaufort Sea and in the north of Ellesmere Island. Areas between the East Siberian Sea and the North Pole are mainly affected by transpolar drift. The Beaufort Sea and the north of Ellesmere Island are within the Beaufort Gyre. Therefore, the possible reason may be that the large Euclidean distance was due to the long ice displacement carried by transpolar drift. Moreover, the ices rotate in the Beaufort Gyre zone, which leads to high angle errors of retrieved ice motions and large cosine distances.

Seasonal and Spatial Variations of The Product Accuracy
As can be seen from Table 3 and Figure 3, the accuracies of Ifremer_SSMIS_H, Ifre-mer_SSMIS_V, and Ifremer_ASCAT were similar to that of Ifremer_Merged. Therefore, these three single-sensor Ifremer products are not included in the present seasonal and spatial variation analysis. Table 6 shows the seasonal errors of the other six annual products. It was found that the freezing season accuracy of most products retrieved by microwave data was significantly better than their melting season accuracy. This may be because the melting ice has a significant impact on the acquisition of bright temperature data and backscattering data, which makes sea ice tracking more difficult. OSI-407-a retrieved by AVHRR performed better in the melting season than the freezing season. This may be because optical remote sensing images are less affected by melting ice. Furthermore, the visible and infrared bands of AVHRR were both used in the melting season, while only the infrared band of AVHRR was used for the freezing season. 5. Furthermore, from Figure 10, it may be seen that the Euclidean distances of almost all products between the East Siberian Sea and the North Pole were larger than those near the Beaufort Sea and in the north of Ellesmere Island. However, from Figure 11, it may be seen that the cosine distances located between the East Siberian Sea and the North Pole were smaller than those near the Beaufort Sea and in the north of Ellesmere Island. Areas between the East Siberian Sea and the North Pole are mainly affected by transpolar drift. The Beaufort Sea and the north of Ellesmere Island are within the Beaufort Gyre. Therefore, the possible reason may be that the large Euclidean distance was due to the long ice displacement carried by transpolar drift. Moreover, the ices rotate in the Beaufort Gyre zone, which leads to high angle errors of retrieved ice motions and large cosine distances.

Seasonal and Spatial Variations of The Product Accuracy
As can be seen from Table 3 and Figure 3, the accuracies of Ifremer_SSMIS_H, Ifremer_SSMIS_V, and Ifremer_ASCAT were similar to that of Ifremer_Merged. Therefore, these three single-sensor Ifremer products are not included in the present seasonal and spatial variation analysis. Table 6 shows the seasonal errors of the other six annual

Impact of Data Source, Retrieval Algorithm, and Time Interval on Accuracy
In order to compare the products with different data sources, two compared groups of eight products derived from the same algorithms with the same time intervals were selected. Table 9 shows the accuracy of the two compared groups. In the OSI SAF compared groups, OSI-405-c-merged in freezing season were compared to keep the consistency in the data period. It can be found from the two groups that products derived from AMSR2 were better than those derived from ASCAT when using either algorithm. The results of the Ifremer group show that the order of ice drift errors is AMSR2 < ASCAT < SSMIS_V < SSMIS_H < merged. This indicates that the merging algorithm can increase spatial coverage, albeit with a decrease in accuracy. Table 6. Seasonal errors of six annual products.

Product
The Freezing Season The Melting Season   Among all eleven products, OSI-405-c_AMSR2 and Ifremer_AMSR2 (freezing season) apply different algorithms but use the same time intervals, data sources, and data periods.
From Tables 5 and 9, OSI-405-c_AMSR2 has a speed MAE of 1.01 km/d and an angle MAE of 11.91 • , while Ifremer_AMSR2 has a speed MAE of 1.15 km/d and an angle MAE of 17.29 • . This shows that the accuracy of OSI-405-c_AMSR2 is higher than that of Ifremer_AMSR2 in terms of speed and angle MAEs, especially in angle. This result demonstrates that CMCC outperforms MCC. The reason is that MCC adopted by Ifremer_AMSR2 is a block-based drift retrieval method that searches the best correlation between two sub images at the pixel level. CMCC adopted by OSI-405-c_AMSR2 searches the candidate vector in a continuous region of the image plane, which can derive the ice drift at the subpixel scale and reduce quantization noise.
In order to fully evaluate the impact of time interval on product accuracy, products with different temporal resolutions derived by the same algorithms and the same data sources should be compared. Ifremer_Merged has products with 3 d and 6 d time intervals, and Ifremer_AMSR2 has 2 d, 3 d, and 6 d products. However, all daily products have no corresponding products with multiday time lags. Therefore, in this study, the two-day products were calculated based on the corresponding daily products by using the trajectory reconstruction method described in Section 2.3.2. In the three daily products evaluated in this study, NSIDC_Pathfinder and NSIDC_AMSR2 were used. OSI-407-a was not used as it has poor spatial continuity and thus cannot be employed in trajectory reconstruction. The estimation results of Ifremer_AMSR2, Ifremer_Merged, NSIDC_Pathfinder, and NSIDC_AMSR2 with different time intervals are shown in Table 10. It is obvious that the error of long-time resolution products is lower than that of short-time resolution products, regardless of angle or speed. The reason may be that the impact of noise in images on derived motions decreases with the increase of time interval. Furthermore, by comparing Tables 3 and 10, it is found that the speed MAEs of NSIDC_Pathfinder (two-day) and NSIDC_AMSR2 (two-day) are 1.44 km/d and 2.15 km/d, still higher than that of OSI-405-c_Merged (i.e., 1.38 km/d) and Ifremer_AMSR2 (i.e., 1.15 km/d) with two-day lags. Moreover, the angle MAE of two-day NSIDC_AMSR2 (i.e., 18.71 • ) is higher than that of two-day OSI-405-c_Merged (i.e., 15.07 • ) and two-day Ifremer_AMSR2 (i.e., 17.09 • ). The possible reason may be that NSIDC_Pathfinder does not include AMSR2 as its data source and NSIDC_AMSR2 only use one channel (i.e., 36.5GHz V pol.) of AMSR2 to retrieve sea ice motion vectors.

Discussion
The above results show that the absolute error of the ice drift speed is positively related with drift speed. Given that sea ice drift is greatly affected by wind forcing [39], this study further analyzed the relationship between the speed MAE of ice drift and the wind speed. OSI-405-c_Merged, NSIDC_Pathfinder, and Ifremer_AMSR2 were selected as the sea ice drift datasets because they performed better than other products based on the above results. Hourly 10 m wind data of ERA5 were used as the wind speed dataset to generate the wind speed information with the same temporal resolution as that used in sea ice drift datasets. The ERA5 is a widely used reanalysis dataset, produced by the European Centre for Medium-Range Weather Forecasts (ECWMF) [36,40]. Figure 12 shows scatter plots between the ERA5 10 m wind speed and absolute error of sea ice speed from OSI-405-c_Merged, NSIDC_Pathfinder, and Ifremer_AMSR2. The figure shows a strong positive correlation between wind speed and the speed error of ice drift, as expected, which indicates the high reliability of the three sea ice drift products. Furthermore, it was found that with an increase of wind speed, NSIDC_Pathfinder showed the largest error increase. The reason for this may be that NSIDC_Pathfinder employs NCEP/NCAR wind forecasts to refine satellite-derived sea ice motions, which means that its ice drift speed and speed errors are strongly related to wind speed.
speed. OSI-405-c_Merged, NSIDC_Pathfinder, and Ifremer_AMSR2 were selected as the sea ice drift datasets because they performed better than other products based on the above results. Hourly 10 m wind data of ERA5 were used as the wind speed dataset to generate the wind speed information with the same temporal resolution as that used in sea ice drift datasets. The ERA5 is a widely used reanalysis dataset, produced by the European Centre for Medium-Range Weather Forecasts (ECWMF) [36,40]. Figure 12 shows scatter plots between the ERA5 10 m wind speed and absolute error of sea ice speed from OSI-405-c_Merged, NSIDC_Pathfinder, and Ifremer_AMSR2. The figure shows a strong positive correlation between wind speed and the speed error of ice drift, as expected, which indicates the high reliability of the three sea ice drift products. Furthermore, it was found that with an increase of wind speed, NSIDC_Pathfinder showed the largest error increase. The reason for this may be that NSIDC_Pathfinder employs NCEP/NCAR wind forecasts to refine satellite-derived sea ice motions, which means that its ice drift speed and speed errors are strongly related to wind speed. According to the research on the general user needs for satellite sea ice data, more than 91% of users indicated that satellite-derived sea ice drift information is useful for their research [41,42]. The main focuses of their research are model validation and data assimilation [41,42]. As model validation and simple data assimilation (e.g., nudging) generally require high accuracy of the sea ice drift information [28], OSI-405-c_Merged and Ifremer_AMSR2 are recommended for such applications because the overall accuracies of these products are high. For advanced data assimilation techniques which employ errors to adjust the weighting coefficients, all the products evaluated in this study could be used as long as the errors of the products are appropriately taken into account [28]. However, NSIDC Pathfinder, with large spatial coverage, is preferable in such cases. Regarding the requirements of satellite-derived sea ice drift products, the World Meteorological Organization states that the target accuracy is 1km/d on a weekly basis, and the spatial resolution should be 5 km [43,44]. From the above results, currently, the studied products achieved accuracies of 1km/d on a weekly basis. However, there are still many ways to further improve the accuracy of sea ice motions, such as by merging vectors retrieved based on the CMCC algorithm using 89 GHz channels of AMSR2. Moreover, there is still a great demand for improving the spatial resolution of satellite-derived sea ice drift products. One possible solution is to develop more applicable methods for high resolution SAR and optical datasets.

Conclusions
In this paper, eleven satellite-derived sea ice motion products, namely, four from OSI SAF, two from NSIDC, and five from Ifremer, have been evaluated by using 647 buoys with high positioning accuracy in the Arctic from 2018 to 2020. The conclusions are as follows: (1) Among the eleven products, NSIDC_Pathfinder, Ifremer_AMSR2, and OSI-405-c_Merged performed well. NSIDC_Pathfinder had the highest angle accuracy among all products but its speed MAE was large. Compared with other products, it had higher temporal (i.e., daily) and spatial resolution (i.e., 25 km) and wider spatial coverage. Ifremer_ASMR2 had the highest speed accuracy, but its angle accuracy and trajectory reconstruction results performed relatively poorly. The spatial resolution of this product is 31.25 km. OSI-405-c_Merged showed the best ability in trajectory reconstruction. Its angle MAE is about 0.1 • larger than that of NSIDC_Pathfinder and its speed MAE is 0.2 km/d greater than that of Ifremer_AMSR2. The spatial resolution of this product is 62.5 km, i.e., the lowest of these three products. (2) In terms of seasonal and spatial variations of the product accuracy, the accuracy of the freezing season was significantly better than that of the melting season for most products derived by microwave sensors. However, the freezing season accuracy of the optical-sensor product (i.e., OSI-407-a) was slightly lower than its melting season accuracy. In addition, the speed MAEs in regions where ice moves faster (i.e., KBS and EG) were greater. The angle MAEs of CBS and EG were higher as ice movements are complex in these regions. Overall, most products performed worst in EG.
(3) Product accuracy can be affected by the data sources, extraction algorithms, and time intervals. The accuracy achieved by different data sources, from good to bad, may be ordered as follows: AMSR2 > ASCAT > SSMIS_V > SSMIS_H. Merging from singlesensor products may not improve accuracy. Moreover, CMCC is superior to MCC in speed and angle accuracy, especially in angle. Furthermore, the accuracies of the products with long-time resolution are better than those with short-time resolution.
Overall, the results of this study indicate that all products are of practical value for model verification and data assimilation if uncertainties and errors are given in an appropriate way. In addition, from the analyses on how product accuracy is affected by season, region, data source, retrieval algorithm, and time interval, it was found that there are many possibilities for the enhancement of sea ice drift monitoring based on remotely sensed data. This paper provides some references for future improvements of these and similar products.