Southern Ocean sea ice concentration budgets of five ocean-sea ice reanalyses

In this study, sea ice concentration (SIC) budgets were calculated for five ocean-sea ice reanalyses (CFSR, C-GLORSv7, GLORYS12v1, NEMO-EnKF and ORAS5), in the Southern Ocean and compared with observations. Benefiting from the assimilation of SIC, the reanalysis products display a realistic representation of sea ice extent as well as sea ice area. However, when applying the SIC budget diagnostics to decompose the changes in SIC into contributions from advection, divergence, thermodynamics, deformation and data assimilation, we find that both atmospheric and oceanic forcings and model configurations are significant contributors on the budget differences. For the CFSR, the primary source of deviation compared to other reanalyses is the stronger northward component of ice velocity, which results in stronger sea ice advection and divergence. Anomalous surface currents in the CFSR are proposed to be the main cause of the ice velocity anomaly. Furthermore, twice the mean ice thickness in the CFSR compared to other reanalyses makes it more susceptible to wind and oceanic stresses under Coriolis forces, exacerbating the northward drift of sea ice. The C-GLORSv7, GLORYS12v1 and NEMO-EnKF have some underestimation of the contribution of advection and divergence to changes in SIC in autumn, winter and spring compared to observations, but are more reasonable in summer. ORAS5, although using the same coupled model and atmospheric forcing as C-GLORSv7 and GLORYS12v1, has a more significant underestimation of advection and divergence to changes in SIC compared to these two reanalyses. The results of the SIC budgets of five ocean-sea ice reanalyses in the Southern Ocean suggest that future reanalyses should focus on improving the modelling of sea ice velocities, for example through assimilation of sea ice drift observations.


Introduction
Antarctic sea ice plays a key role in the global climate system (Thomas and Dieckmann 2010). However, sea ice studies in Southern Ocean are limited by both the availability of observations and the reliability of numerical models. On the one hand, in-situ observations of sea ice are very scarce due to harsh climatic conditions. Therefore, satellite observation is the primary reliable data source to obtain sea ice concentration (SIC). For other parameters such as sea ice thickness and sea ice velocity, the reliability of satellite observations remains challenging (e.g., Giles et al. 2008;Xie et al. 2011;Kurtz and Markus 2012;Huntemann et al. 2014). On the other hand, the results of current climate models in Antarctica deviated from the fact of overall increase of the Antarctic sea ice extent (SIE) observed over the past four decades (Turner et al. 2013;Zunz et al. 2013;Mahlstein et al. 2013), even though it experienced a dramatic decline over the period 2015 to 2017 (Parkinson 2019;Eayrs et al. 2021).
As a reconstruction of past ocean and sea ice states, ocean reanalysis, including sea ice parameters (referred to as iceocean reanalysis), provides a valuable database for studies of polar sea ice. In comparison to coupled climate models, ice-ocean reanalyses integrate ocean and sea ice observations through data assimilation and indirectly integrate meteorological observations through fixed atmospheric forcing coming most often from atmospheric reanalyses (Dee et al. 2014;. Compared with in-situ or satellite observations, ice-ocean reanalysis is not only continuous in time and space, but also shows inherent ocean model physics (Storto et al. 2019). In view of the above advantages, ice-ocean reanalysis is widely used to study sea ice trends and sea ice-atmosphere interactions in polar regions (e.g., Ponsoni et al. 2018;Spreen et al. 2011;Jaiser et al. 2016), or as initial and boundary conditions for sea ice and atmosphere forecasting models (e.g., Guemas et al. 2016), and numerical weather prediction models (e.g., Dee et al. 2011). However, due to differences in model formulation, boundary conditions and data assimilation approaches, there is also significant spread between reanalysis products themselves . Therefore, a systematic inter-comparison of reanalysis products is essential to better understand their advantages as well as limitations.
Assessments of Southern Ocean ice-ocean reanalysis products are still relatively scarce. Uotila et al. (2019) conducted a detailed comparison of nine reanalysis products in the Southern Ocean, focusing on sea ice edge, SIC and sea ice thickness. They concluded that the SIE and sea ice area (SIA) of the reanalyses were in good agreement with observations, except for one product, SODA, which was a clear outlier, while the sea ice thickness of the reanalyses was more spread. Shi et al. (2021) evaluated sea ice thickness for four reanalyses in the Weddell Sea based on observations from altimeters and in-situ measurements and revealed that the accuracy of these reanalyses is inconsistent with reference to observations from altimeters and in-situ measurements. Both studies are valuable in broadening our understanding of reanalysis in the Southern Ocean, but the complexity of the models makes it difficult to identify the links between errors in the variables and the sources of model bias from the analysis of individual variables (Lecomte et al. 2016). It is therefore expected that further diagnostics of the ice-ocean reanalyses will identify whether the model bias originates from the misrepresentation of dynamic or thermodynamic processes and will guide research to reduce these biases accordingly. Otherwise, the assimilation of SIC may lead to even poorer estimates of sea ice state if no adjustment is made (Dulière and Fichefet 2007), which is clearly contrary to the original purpose of constructing the ice-ocean reanalysis product. Holland and Kwok (2012) proposed a method for calculating the SIC budget in observations by decomposing the change in SIC into advection, divergence and residual terms, with the residual term containing the thermodynamic freezing/melting and mechanical redistribution effects. This approach has recently been used to diagnose climate models (Uotila et al. 2014;Lecomte et al. 2016;Holmes et al. 2019) and to investigate the effects of atmospheric forcing on sea ice simulation results (Barthélemy et al. 2018). In addition, a budgeting method for sea ice thickness/volume has been proposed (Holland et al. 2010), but due to the high uncertainties in sea ice thickness observations in the Southern Ocean, the application of this method is currently limited to comparisons between coupled climate models, such as the ones participating in the Coupled Model Intercomparison Project Phase 5 (Schroeter et al. 2018) and Phase 6 .
The objectives of this study are (1) to investigate whether the dynamic and thermodynamic contributions to SIC changes in Southern Ocean presented by the ice-ocean reanalyses are realistic or not. (2) What are the sources of deviation in terms of physical processes and forcings if they are deviated from the observations? Through this diagnosis, we expect to provide an improved understanding of the Southern Ocean ice-ocean reanalysis products, and whether they simulate realistic SIC patterns for the right reasons. The availability of observed SIC budgets (Holland and Kimura 2016) now makes this objective feasible.

Ice-ocean reanalyses
Five ice-ocean reanalyses selected for evaluation in this study are summarised in Table 1. These reanalyses were chosen because they output SIC and sea ice velocity at the daily frequency. The National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) is a global, high resolution reanalysis product (Saha et al. 2010). The CFSR version 1 covers the period from 1979 to 2010 and extends to the present in its second version (CFSv2). Both versions use the same ocean model, i.e., the Modular Ocean Model version 4p0d (MOM4) coupled with the Sea Ice Simulator from Geophysical Fluid Dynamics Laboratory (GFDL SIS) (Saha et al. 2014). It is worth noting that although the CFSR uses a coupled atmosphere-ocean model, it is in fact an uncoupled reanalysis because the interaction of the atmosphere and ocean is not directly used (Alkama et al. 2020). Sea ice thickness in each grid cell is divided into five categories, i.e., 0-0.1, 0.1-0.3, 0.3-0.7, 0.7-1.1 m, and thicker than 1.1 m.
The CMCC Global Ocean Physical Reanalysis System (C-GLORS) is a high-resolution and eddy-permitting reanalysis  and its latest version was used in our study (C-GLORSv7). C-GLORSv7 is based on the Nucleus for European Modelling of the Ocean (NEMO) ocean model, coupled with the Louvain-la-Neuve Sea Ice Model version 2 (LIM2) sea ice model, assimilating satellite-derived SIC through a 6-h nudging scheme in the Southern Ocean and implementing a large-scale bias correction scheme to avoid ocean model biases. Sea ice thickness in C-GLORSv7 is weakly constrained by a 15-day nudging scheme in the Arctic Ocean but not in the Southern Ocean, which is mainly due to the lack of sea ice thickness observations or validated reanalysis in the Southern Ocean.
The Copernicus Marine Environment Monitoring Service (CMEMS) global ocean 1°/12° reanalysis product (GLO-RYS12v1), which is eddy-permitting, has the highest spatial resolution among these reanalyses. The SIC is assimilated by using the singular extended evolutive Kalman filter (Brasseur and Verron 2006). Similarly to C-GLORSv7, a 3D-VAR large-scale bias correction is performed in GLORYS12v1. With the development of the reanalysis system, the most recent version is a significant improvement over previous ones, particularly in terms of global heat and salt contents (Lellouche et al. 2021).
The NEMO-EnKF published by Massonnet et al. (2013) was the first SIC-constrained reanalysis of Antarctic sea ice thickness over the last few decades, and was therefore included in our assessment as an important benchmark. The SIC is assimilated into the NEMO-LIM2 model by means of ensemble Kalman filtering in the reanalysis, but sea surface temperature is not, which is different from that of other four reanalyses. A "FREE" run without assimilation was also carried out simultaneously and the positive effect of assimilating SIC on improving sea ice extent and thickness was verified by comparing these two experiment results to independent ship-based estimates of sea ice thickness. Another difference between NEMO-EnKF and other reanalysis products is that the atmospheric forcing used is NCEP/NCAR (Kalnay et al. 1996), which is reported to be questionable (Bromwich et al. 2007;Massonnet et al. 2013) mainly due to a reported positive bias in wind speed over the Southern Ocean (Li et al. 2013). The ECMWF Ocean ReAnalysis System 5 (ORAS5; Zuo et al. 2019), containing five ensemble members generated by the perturbation of initial conditions, observations and forcings, is another high-resolution and eddy-permitting reanalysis analysed in this study. The SIC in ORAS5 is almost the same as that of ORAP5 (Zuo et al. 2015), which is the prototype of ORAS5 and has been proved to be a relatively good representative of Antarctic sea ice (Uotila et al. 2019). Similarly to the other reanalyses, we used the unperturbed member of ORAS5, since the unperturbed member is usually the one disseminated and validated and used for comparison (Zuo et al. 2019). There is almost no difference in the SIC budget between the members (not shown). The SIC assimilation approach used by ORAS5 is the First-Guess at Appropriate Time (FGAT), a three-dimensional variational assimilation (3D-Var) method.
Comparing the models and coupling schemes used in these reanalyses, except for CFSR, the remaining four all use NEMO-LIM2 model and have a single sea ice thickness category. Three of them (C-GLORSv7, GLORYS12v1 and ORAS5) are developed in the framework of the MyOcean and MyOcean2 projects and have many commonalities, e.g., they are all high-resolution, eddy-permitting and using ERA-Interim (Dee et al. 2011) as atmospheric forcing. In the Sourthern Ocean, all these five reanalyses assimilated only SIC, the other two key sea ice variables, sea ice drift and sea ice thickness, were not assimilated.

Observations from satellites
Two satellite-derived sea ice velocity products are used in the SIC budget calculation (Table 2). One is derived from 36-GHz channel brightness temperature data of Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E), with an average spatial resolution of 60 km (Kimura et al. 2013), and the other is Polar Pathfinder Daily 25 km EASE-Grid Version 4 product from the National Snow and Ice Data Centre (NSIDC; Tschudi et al. 2019). Both of them are achieved with maximum cross correlation techniques, the difference is that the latter one derived from various sources is merged to produce daily sea ice motion. The initial version of the NSIDC product was reported to have significantly underestimated ice speeds compared to buoy observations (about 40%; Heil et al. 2001), while Kimura et al. (2013) has an overestimation of 5%. Although the algorithm in Version 4 has been significantly improved, with an overall increase in sea ice speed of 0.5-2 cm/s compared to Version 3 (Tschudi et al. 2019), we still corrected the NSIDC ice drift velocity by multiplying the meridional speed with a constant factor equals to 1.357 (Haumann et al. 2016).
Daily SIC derived from AMSR-E (Cavalieri et al. 2014) brightness temperature data by using Enhanced NASA Team (NT2) algorithm is adopted in the calculation of observational SIC budgets, as that in Holland and Kimura (2016). Four satellite-derived SIC observations are also used as references when comparing SIC between the AMSR-E and reanalyses (Table 2)  (ASI) algorithm. The monthly climatologies of these four observation products are calculated from their daily outputs. All data described above, including ice-ocean reanalyses and satellite observations, are regridded onto the Kimura et al. (2013) ice drift product's 60 km common grid by linear interpolation for the comparability. Their evaluation is limited to a common period (from February 2003 to October 2009) due to the availability of the NEMO-EnKF reanalysis and Kimura et al. (2013) ice drift data.

SIC budget
According to Holland and Kwok (2012), the local temporal rate of change of SIC ( A ) is governed by dynamic and thermodynamic processes. However, for reanalysis products that have assimilated SIC, the contribution of the data assimilation to SIC changes cannot be ignored. Thus, the SIC budget equation can be expressed as: where is the sea ice velocity. The right-hand side of the equation represents all non-resolved changes in SIC, including the changes from freezing/melting ( f ), mechanical ice redistribution ( r ), such as ridging and rafting, and data assimilation ( s ). After integrating and reorganizing, Eq. (1) can be rewritten as: i.e. the net change in SIC from the beginning to end of a time period consists of the integration of advection, divergence, thermodynamic, deformation and data assimilation components over the same time period. For the sake of simplicity, the terms in Eq. (2) are denoted as change or dadt , advection or adv , divergence or div and residual or res from left to right. Apparently, the positive values of these terms correspond to an increase in SIC, and for observations the data assimilation term is always zero. Following Holland and Kwok (2012) and Haumann et al. (2016), a 3 × 3 cell filter is used to smooth out grid-scale noise in satellite-derived ice drift field. For consistency, this low pass filter is also applied in reanalyses, although the noise mentioned above does not exist inside model output (Uotila et al. 2014).
The daily and seasonal calculation of each term adopts the same approach as Holland and Kimura (2016). Specifically, daily changes are calculated by central differencing of the SIC fields for the day before and after; the advection and divergence are firstly calculated on daily basis and then averaged over the same 3-day periods to be consistent with the daily SIC change; then, the residual term f − r + s is obtained by subtracting adv and div from dadt.

Spatial and seasonal divisions
For easier understanding of how SIC budgets and sea ice velocity in reanalyses differ from observations regionally, the Southern Ocean is zonally divided into five sectors between longitudes 20° E, 90° E, 160° E, 120° W and 50° W (Fig. 1). Moreover, the dynamics of sea ice in areas close to the coast, driven by easterly winds and with high SIC, is very different from the one in areas further offshore, driven by westerly winds and with relatively low SIC. The study area, bounded by 40°S, is therefore also divided into "interior" and "exterior" areas as in Lecomte et al. (2016). The border of these two areas is based on the February 2003 to October 2009 mean ERA-Interim 10 m wind, following the change in wind direction to better distinguish between the westerly and easterly wind regimes. Seasonal divisions in this study are defined in Fig. 2.

Results and discussion
As the most basic diagnostic, the Sect. 3.1 compares SIC in the reanalysis products with the observed SIC. The results of the SIC budget are presented in the Sect. 3.2. Sea ice velocity and sea ice thickness from the reanalyses are briefly compared with observations in Sects. 3.3 and 3.4, respectively, to further explore the reasons why the SIC budgets in the reanalyses differ from the observed. Figure 2 displays the mean seasonal cycle of some key measures related to SIC. Not surprisingly, thanks to the SIC data assimilation, the sea ice extent (defined as the integral of grid cells areas where the SIC is larger than 15%) and area (the integral of grid cells areas multiplied by the SIC in each grid cell) in sea ice reanalyses are generally consistent with that observed. Figure 2 also shows the root mean square error (RMSE) between the reanalyses and the four satellite observations used for assimilation, defined as where x obs i denotes the observation at month i.

Sea ice concentration
The sea ice extent of AMSR-E agrees well with the other observations (RMSE = 3.8 × 10 5 km 2 ), but its marginal ice zone (MIZ; area covered with SIC between 15 and 90%) area is too small and its pack-ice area (area covered with SIC higher than 90%) is too large, which is responsible for its larger total sea ice area than the other observations. This is attributed to the retrieval algorithm used by AMSR-E which overestimates the percentage of sea ice within a pixel cell (Beitsch et al. 2012). The sea ice area of CFSR is quite close to that of AMSR-E, the difference is that CFSR has a more realistic MIZ and pack-ice area. C-GLORSv7 has lower than other datasets' sea ice extent and area in summer, but this does not occur in other seasons, which results in a higher February to March sea ice extent growth rates (Fig. 2a). GLORYS12v1 has the closest sea ice area (RMSE = 6.4 × 10 5 km 2 ) to observations of all reanalyses, with all months within the observed range except from January to March, when sea ice extent is slightly below observations. NEMO-EnKF overestimates the sea ice extent in autumn and winter but matches in the other months and has a similar high MIZ/pack-ice ratio increase rate in spring as GLORYS12v1. Finally, ORAS5 is very close to C-GLORSv7, except for a slightly higher pack-ice area and sea ice area in December.
As shown in Fig. 2, the sea ice extent and area are relatively consistent across all datasets but diverge more on the MIZ and pack-ice area, even among the four observation datasets. As mentioned above, this may stem from systematic bias due to differences in retrieval algorithms. However, the MIZ and pack-ice areas of these reanalyses are almost always within the observed range, in slight contrast to Uotila et al. (2019), who showed in their study that sea ice data assimilation had little effect on the MIZ/pack-ice ratio. This slightly different result may be due, on the one hand, to the different time frame of the study, and on the other hand, to the fact that the reanalyses used in Uotila et al. (2019) have been updated.
Although AMSR-E has a too low MIZ/pack-ice ratio, which is potentially related to the NT2 retrieval algorithm it used (Beitsch et al. 2012), it is still used to build the observation-based SIC budget in this study, due to the supplementary experiments of Holland and Kimura (2016), which revealed that the SIC budget is insensitive to the uncertainty in the SIC dataset.  Figures 3,4,5,6 show the SIC budget of five reanalyses compared to observations in different seasons, the two observed SIC budgets are denoted as OBS-K and OBS-T, being calculated from Kimura et al. (2013) and Tschudi et al. (2019) sea ice drift data, respectively. The Student's t-test was applied to test whether the time series of the SIC budget component were significantly different between the reanalyses and observations. From observation-based results, it appears that in all seasons, advection processes contribute to local increase in SIC at the ice edge, where high air temperatures are causing ice melting. The divergence causes a reduction in SIC, resulting in more open water and favourable conditions for sea ice growth.

SIC budgets in reanalyses
To quantify the contribution of each component to the SIC change, Table 3 provides statistics on the Antarctic-wide area integrals of the SIC budget components in four seasons. It should be noted that sea ice changes in the reanalyses are larger than those in the two observational datasets ( dadt columns in Table 3), because dadt is only integrated over the grid points where the ice drift data are available, and more ice drift data equals zero in the two observations than from the five reanalyses. Since dynamical processes only spatially redistribute the existing SIC, the contribution of dynamical processes to the change in total sea ice area should be zero. However, due to the inaccuracy and missing sea ice velocity data in the observational dataset, there is an imbalance in the observed advection and divergence area integrals (Holland and Kimura 2016).
During the autumn, sea ice is growing along the Antarctic coast, particularly in the Ross and Weddell Seas, except for the western Weddell Sea (the dadt column in Fig. 3), where a high SIC is maintained all year round. All five reanalyses reproduce the spatial patterns of sea ice growth well and are in good agreement with the observations. In contrast, there are more differences when the SIC change is decomposed into the advection, the divergence and the residual Grey crosses mark statistically significant differences between the reanalyses and both two observational products at the 95% confidence level according to the t-test components. In the reanalyses, the advection, in addition to causing an increase in SIC at the outer edge region as in the observations, also causes a decrease in SIC along the coast (the adv column in Fig. 3). As mentioned in Holland and Kwok (2012) and Uotila et al. (2014), observationbased SIC budgets are unable to resolve sea ice drifts in sub-pixel geometries near the coastline due to low resolution. In addition, the applied smoothing plays considerably levels the small-scale coastal features out. Therefore, the difference between reanalyses and observations in terms of coastal advections cannot be evaluated at this point. Apart from CFSR, all reanalyses show a lower than observed advection-induced sea ice increase as a proportion of the total increase in SIC (Table 3).
The divergence is the term that features the most differences between reanalyses and observations (the div column in Fig. 3). All reanalyses are significantly different from observations in almost all regions, with CFSR having overall too strong divergence, while the other reanalyses have too weak divergences. Specifically, CFSR overestimates the reduction in SIC due to divergence in the western Weddell Sea, the Amundsen-Bellingshausen Seas, the Ross Sea, the Dumont d'Urville Sea and the Prydz Bay. This results in the CFSR producing excessive thermodynamic sea ice growth at these areas to offset the sea ice reduction due to divergence (see the residual of the CFSR in Fig. 3). The other four reanalyses all underestimate the divergence, both in terms of magnitude and the extent, especially in the Ross Sea, the Weddell Sea and the Prydz Bay. However, the overestimation of advection and divergence in the CFSR and the underestimation in the other reanalyses are balanced and eventually result in a reasonable residual term.
On the average, sea ice in the Southern Ocean grows at lower latitudes in winter than in autumn (the dadt column in Fig. 4). Again, the reanalyses reproduce the overall tendency of SIC but fail to accurately capture the contribution of advection and divergence to SIC change. The CFSR transports too much sea ice to north via advection, resulting in more sea ice melting at the ice edge to maintain a reasonable SIC distribution, while excessive divergence leads to a reduction in SIC in the inner ice pack, resulting in excessive sea ice growth (residual of the CFSR in Fig. 4). For the remaining reanalyses, although the advection looks relatively accurate, it does not agree well with the observations in terms of sea ice fluxes in the inner ice zone, as the t-test shows them to be statistically significantly from different distributions. Despite the large discrepancies in the divergence terms based on the two different sea ice velocity products, they both indicate strong divergence in the Ross Sea, the Weddell Sea, the Dumont d'Urville Sea and the Prydz Bay. However, C-GLORSv7, NEMO-EnKF and ORAS5 underestimate the divergence in these regions, GLORYS12v1 has close to observed divergence in the the Dumont d'Urville Sea and the Prydz Bay but not in the Ross Sea and the Weddell Sea. This is also evidenced in Table 3, where in each reanalysis the proportion of SIC change accounted for by divergence is lower than observed.
The NEMO-EnKF product shows high sea ice convergence (i.e. positive value of div term) in the West Indian sector, the East Indian sector and the Amundsen-Bellingshausen Seas unlike other systems ( div column in Fig. 4). In addition, the convergence in the West Indian sector from NEMO-EnKF is balanced by a reduction in SIC due to advection, resulted in a relatively consistent with observations' residual in this sector. For ORAS5, the contribution of advection and divergence to sea ice change is very small, at only 5% and -2% respectively, a proportion essentially the same as its autumnal counterpart.
In spring, sea ice starts to melt from the edge of the ice zone at lower latitudes (Fig. 5). For both observations and reanalyses, except NEMO-EnKF and ORAS5, the overall spatial distributions of sea ice advection as well as divergence are similar to those of winter, and the area integrals are also close to the winter values (Table 3). For the CFSR, the contribution of advection and divergence to SIC change is still anomalously excessive compared to observations. Additionally, the convergence of sea ice in the eastern shore of the Bellingshausen Sea does not occur in the CFSR but it does in the other four reanalyses ( div column in Fig. 5), which is possibly due to a bias in the direction of sea ice velocity in the CFSR (Fig. 7). The spatial pattern as well as the area integration of the SIC budget components of C-GLORYSv7 and NEMO-EnKF are quite consistent during this season. GLORYS12v1 is closest to OBS-T in terms of the proportion of advection and divergence contributions to SIC change. The performance of advection and divergence in ORAS5 is better compared to autumn and winter as its percentage has improved (Table 3), but is still significantly underestimated compared to observations. During the summer months, most of the sea ice melts as the temperature rises. The advection in the CFSR transports large amounts of sea ice to north, where it melts, being  (100) compensated by the divergence causing an excessive reduction of sea ice (Fig. 6). In addition, the CFSR produces excess sea ice growth to compensate for the divergencedriven sea ice reduction in the western Weddell Sea. The other four reanalyses show a dynamical contribution to sea ice change very close to that observed (Table 3), which can be explained by the fact that during summer sea ice change is mainly thermodynamic melting, with a small contribution of about 10% from dynamics. As shown in Fig. 6, the grids where advection and divergence are statistically significantly different from observations in these four reanalyses are much fewer compared to the other seasons. Additionally, Fig. 7 Seasonal mean Antarctic sea ice velocity vectors (arrows) and speed (color shades) for five reanalyses (5 bottom rows) and two observational products (2 upper rows), over the study time period the patterns of the residual terms are also basically similar to observations, for example, all showing areas where SIC keeps increasing during the summer months, such as in the western Weddell Sea, the eastern coast of the Ross Sea, and the western Dumont d'Urville Sea (the res row in Fig. 6).

The contribution of ice velocity to the deviations of SIC budget
The previous diagnostic results for the SIC budget have shown that the difference between reanalyses and observations SICs lies mainly in the dynamically driven sea ice variability, thus, Fig. 7 provides a further insight into sea ice velocity. Comparison of the two observations shows that although the Tschudi et al. (2019) ice velocity is corrected by by enlarging the meridional component, its velocity magnitude remains clearly smaller than that of the Kimura et al. (2013). Sea ice drift speeds of all reanalyses are underestimated in the marginal ice areas of the Ross and Weddell Seas in summer and autumn, and slightly underestimated in the Amundsen-Bellingshausen Sea in winter. The five reanalyses all exhibit relatively large sea ice speeds along the coasts of the Weddell Sea, the West Indian sector, and the East Indian sector, which could be attributed to the high model resolution of the reanalyses (excepte NEMO-EnKF) and thus the ability to resolve the coastal grid ice drift. Among these reanalyses, it is clear that sea ice velocities of GLORYS12v1 are closest to the observations in the Ross Sea in autumn, winter and spring. Excluding the above differences, Fig. 7 shows that the sea ice drift speeds of the reanalyses are essentially in between the two observational estimates. Therefore, the differences of SIC budgets between the reanalyses and the observations studied here seem to be more sensitive to the direction of sea ice drift. Subsequently, to complement Fig. 7, Fig. 8 displays only the north-south component of sea ice velocity, as the east-west component forms a ring around the Southern Ocean and plays a little role in the annual SIE expansion. Moreover, it can be expected that in general, higher north-south ice drift speeds correspond to the increased SIC at the northern ice edge due to stronger advection, while higher north-south velocity gradients correspond to stronger divergence. In autumn, the predominantly westerly direction of the reanalyzed sea ice drift in the interior area of the Ross Sea and the lack of a northerly component in the exterior area (the first column in Figs.7, 8) result in advection-driven SIC increase in the reanalyses (except GLORYS12v1) mainly in the eastern Ross Sea and in the westernmost part of the East Indian section (Fig. 3). Strong north-south velocity gradients in the observations are responsible for the divergence of this area. The north-south velocity gradient in the CFSR is not prominent compared to observations and other reanalyses in autumn. This also corresponds to the fact that the proportion of divergence-driven sea ice reduction in the CFSR is closest to OBS-K during this season compared to other seasons (Table 3).
In winter and spring, the SIE in the Southern Ocean is larger than in autumn and the sea ice drift direction of the CFSR is more distinctly characterised by a large northward velocity component in the exterior area and a clear northward velocity gradient from interior to exterior area throughout the whole Southern Ocean (Fig. 8). Clearly, this is the reason for the significantly higher-than-observed advectiondriven sea ice increase at the northern edge of the CFSR during these seasons, as well as the excessive divergence. During summer, as sea ice velocities are all low (Fig. 7), the spatial characteristics of the north-south gradient of the northward velocity component do not differ much between reanalyses, with only a slight underestimation in the Ross Sea compared to observations. The northerly sea ice drift in NEMO-EnKF showing higher than other systems' values around the Oates Land in the western Ross Sea, thus leading to the ridging of sea ice in the region ( div of 5). C-GLORSv7, GLORYS12v1 and ORAS5 have pronounced north-south velocity gradients in front of the Ross Ice Shelf, which allows for a higher reduction in SIC by the divergence, and thus resulting in a higher sea ice growth in this region ( div rows in Figs. 3, 4, 5, 6) than NEMO-EnKF.

Differences in sea ice thickness between reanalyses
The SIC, sea ice velocity and sea ice thickness are interdependent to each other and, thus, an analysis of sea ice thickness is useful for further understanding the mass balance of sea ice in the Southern Ocean. Figure 9 shows the spatial distribution of temporally averaged sea ice thickness for each reanalysis and their monthly climatology and the NEMO-EnKF sea ice thickness corrected by ICESat-1 (2003ICESat-1 ( -2008 and ASPeCt (1980ASPeCt ( -2005 observations (Haumann et al. 2016) is considered here as the reference. It is noticeably clear that throughout the Southern Ocean, sea ice in the CFSR is much thicker than in other systems ( Fig. 9a-f), accompanied by a larger decrease of sea ice thickness from the interior to exterior area. The mean sea ice thickness of CFSR reaches its minimum in February, when the NEMO-EnKF is at its maximum (Fig. 9g). The CFSR sea ice thickness continues to increase until November, reaching 1.5 m, while the NEMO-EnKF slowly increases after a rapid decrease in autumn. This difference in mean sea ice thickness and growth of CFSR has also been reported by Huang et al. (2015). C-GLORSv7, GLO-RYS12v1 and ORAS5 have similar spatial distributions of sea ice thickness and they show that the thickest sea ice occurs along the eastern side of the Antarctic Peninsula.
These three reanalyses are closer to the corrected NEMO-EnKF than the original NEMO-EnKF, which overestimates sea ice thickness in the west-central Weddell Sea. Furthermore, similarly to the corrected NEMO-EnKF, GLORYS12v1 and ORAS5 also have slightly thicker sea ice in the interior regions of the West Indian, the East Indian and the Amundsen-Bellingshausen Seas sectors, which results in larger sea ice volumes in these reanalyses than in NEMO-EnKF (Fig. 9h). Although C-GLORSv7 does not constrain the sea ice thickness in the Antarctic, as it does in the Arctic, it has thinner sea ice compared to other reanalyses, which may stem from the adjustments to Fig. 8 Contour maps of seasonal mean northward sea ice drift speed for different products the data assimilation system and uncorrected wind forcing . The differences in sea ice thickness among the five reanalyses are much greater than in SIE and SIA, owing to the fact that none of these models assimilate sea ice thickness and therefore the calculation of sea ice thickness relies mostly on the model. Nevertheless, sea ice thicknesses in the four reanalyses, apart from the CFSR, are generally reasonable.

General discussion
The reanalyses evaluated in this paper assimilated only SIC in the Southern Ocean, with assimilation techniques ranging from nudging to 3D-Var and Kalman filters. As in previous assessments of ice-ocean reanalyses, the assimilation of observed SIC has undoubtedly resulted in more realistic SIE and SIA in the reanalyses and a significant improvement in the variability of SIC compared to the free runs corresponding to the reanalyses Lellouche et al. 2021) and the reanalysis products without assimilated SIC (Uotila et al. 2019). Assimilation of SIC can also reduce biases on the simulation of sea ice thickness by about 20% ). However, errors in simulated sea ice velocity arise mainly from errors in the wind field or ocean currents (Barth et al. 2015), in particular in the Southern Ocean, where the sea ice internal stress is generally small, and the impact of SIC assimilation on sea ice velocity is actually very limited (Stark et al. 2008;Rollenhagen et al. 2009). Therefore, the ensuing discussion of the differences in sea ice velocities between reanalyses and observations mainly focuses on model physics and atmospheric and oceanic forcing. The CFSR stands out in these reanalyses because of the overestimation of advection and divergence on SIC changes, and the preceding results show that its northward ice drift speed and sea ice thickness are also anomalously large compared to other systems. Zhang et al. (2018) examined the 10 m winds as well as surface-ocean current in several iceocean reanalyses, providing evidence for considerable biases in the CFSR sea ice velocity. As they found, the 10 m wind field in the CFSR is in very good agreement with the ERA-Interim reanalysis, however the surface-ocean currents have strong northeastward anomalies between 50° S and 60° S, which basically corresponds to the exterior area in this study. Additionally, the continuous region of statistically significant anomalous northeastward current in the CFSR indicates that this is a systematic error associated with the model. The different behaviour of the 10 m wind field and ocean surface currents in the CFSR can be explained by its non-direct use of atmospheric and oceanic interactions as mentioned before (Sect. 2.1). The excessive northward surface-ocean currents could then be the main reason that drive excessive northward sea ice advection, transporting more sea ice to the northern boundary region and then reduced by melting or data assimilation there. In addition to the adjustment of the SIC by the data assimilation, simultaneously, in order to maintain realistic SIC, new sea ice needs to be continuously frozen to compensate for the reduction in SIC due to divergence. Furthermore, the excessive sea ice thickness in the CFSR causes it to be more susceptible to wind and ocean stresses, resulting in more northward free-drifting sea ice from the predominantly westerly winds and currents in the Antarctic Circumpolar Current. This conclusion can be demonstrated as follows: the equilibrium relationship between the wind stress, ocean stress and Coriolis forces for free-drifting sea ice is (Leppäranta 2011, Chapter 6.1.1) with the symbols explained in Table 4. The general solution of Eq. (3) can be obtaned by solving the following algebraic equations: The wind factor and deviation angle as the function of sea ice thickness then can be obtained by substituting (4) 4 + 2 sin w RNa 3 + R 2 Na 2 2 − Na 4 = 0, the values of the parameters in Table 4 into Eqs. (4)-(5). As shown in Fig. 10, the free drift wind-ice turning angle is quite sensitive to the thickness of the sea ice, and an increase in sea ice thickness from 0.5 to 1.5 m (Fig. 9f) in the Southern Ocean will cause a turning of about 6 degrees more to the north when the wind is from west. Some studies have indicated that the assimilation of SIC may lead to an increase of sea ice thickness (e.g., Tietsche et al. 2013;. Although there are some approaches, for example, to constrain sea ice thickness (e.g., Luo et al. 2021), or to conserve the ice volume (Dulière and Fichefet 2007), that can reduce too thick ice reproduced by the model, such approaches do not seem to have a fundamental effect on the errors caused by external forcing. The comparison of the SIC budgets for the reanalyses other than CFSR shows that, although they all underestimate the contribution of advection and divergence to sea ice change, Figs. 3, 4, 5, 6 demonstrate that overall they are relatively reasonable. It is important to note that in our study, the SIC components of NEMO-EnKF (with NCEP as the atmospheric forcing) and C-GLORSv7 (with ERA-Interim as the atmospheric forcing) were not significantly different. This differs from the findings in Barthélemy et al. (2016), who showed that the model using NCEP clearly have more dynamic contributions to sea ice change than that using DRAKKAR Forcing Set version 5.2, whcih is based on the ERA-Interim. The reason for this difference is mainly that the dynamical processes defined in Barthélemy et al. (2018) include ridging and rafting, and since NCEP has a reported positive bias in wind speed over the Southern Ocean (Li et al. 2013), it is more likly to cause ice ridging.
Our results are therefore not contradictory to them. Additionally, although ORAS5 uses the same atmospheric forcing as C-GLORSv7 and GLORYS12v1, its advection and divergence contributions to sea ice change are most significantly underestimated compared to observations and other reanalyses. This implies that the model configuration also has a considerable influence on the SIC budget. However, more details on model setups currently not available, would be needed for a more in-depth analysis.

Conclusion
In this study we diagnose SIC budgets in five state-of-the-art reanalyses in the Southern Ocean by comparing them with observations, from February 2003 to October 2009. We find that although the reanalyzed SIC fields are in general good agreement with the observed ones, this is the result of error compensations in the various terms governing the concentration budget equation. The results suggest that unrealistically simulated SIC budgets are mainly caused by biases in ice drift speed, as proposed in previous studies (e.g., Uotila et al. 2014;Lecomte et al. 2016;Barthélemy et al. 2018). Furthermore, we emphasise that the meridional component of sea ice drift is the chief cause for biases in the reanalyses.
The CFSR has the largest error in the SIC budget compared to other reanalyses, and we argue that the anomalies in oceanic forcing, i.e. anomalously large northeastward surface-ocean currents, are one of the main causes. However, quantifying the impact of these anomalous currents on sea ice drift and whether they arise from a misrepresentation of atmosphere-ocean interactions (Alkama et al. 2020), or from the ocean model itself, is beyond the scope of this paper and deserves additional in-depth investigation. In addition, how the excessive advection and divergence in the CFSR relates to its approximately twice the observed sea ice thickness is an interesting and critical topic, and we believe it deserves a separate in-depth study. Once it is demonstrated that excessive sea ice thickness in the CFSR is caused by its excessive advection and divergence, a positive feedback of "excessive sea ice drift northwards-excessive divergence-thicker sea ice-sea ice drift more northerly" will be formed, which will be of great implications for the diagnosis of biases in reanalyses and climate models (such as CMIP6).
The C-GLORSv7, GLORYS12v1 and NEMO-EnKF perform relatively well in this study, followed by ORAS5, although they underestimate, to some varying degree, the effects of advection and divergence on the change of the Southern Ocean SIC. However, given that their advection and divergence are generally balanced, and the sea ice thickness is relatively reasonable, we believe that they can be safely recommended to users. In specific, C-GLORSv7 and GLORYS12v1 are reliable in all four seasons, and ORAS5 in (left y-axis) and deviation angle (right y-axis) as a function of sea ice thickness. The negative sign of the angle represents turning left from the wind spring and summer, and all three of them are more reliable in the West and East Indian Ocean sectors than in other sectors. NEMO-EnKF is reliable in seasons except winter, when its excess convergence in the West Indian sector is significantly different from that observed.
Although there is still uncertainty in the observed sea ice velocity data used in this paper, comparing the SIC budgets from the reanalyses with budgets from observations is still useful for providing further diagnostics for ice-ocean reanalysis than for individual variables. Specifically, this diagnosis allows us to better understand what deviations in these variables reveal about the model physics or boundary conditions or model configurations. In addition, the close approximation of actual sea ice extent and area obtained through SIC data assimilation does provide a significant improvement in the hindcast of SIC change ), but does not appear to play a significant role in improving sea ice velocity simulations. Furthermore, we believe that preserving the increments coming from data assimilation to distinguish them from the thermodynamic growth/melt and redistribution can have a positive effect in the process of diagnosing the model.
Biases in the simulation of sea ice velocities can seriously limit the skill of climate models. For example, current climate models are unable to simulate the sea ice expansion in the Southern Ocean, possibly due to systematic biases in sea ice velocities (Sun and Eisenman 2021). Following the assimilation of SIC and the emergence of Antarctic sea ice thickness satellite observations, Luo et al. (2021) assimilated Antarctic sea ice thickness observations into a coupled sea ice-ocean model for the first time and demonstrating the potential to improve the Antarctic sea ice prediction with data assimilation. Given the importance of boundary conditions and the lack of accuracy of winds in the Southern Ocean, another promising approach is to adjust winds by assimilating sea ice velocities (e.g., Barth et al. 2015).
In conclusion, we believe that the next step of improvement should focus on the assimilation of sea ice velocity, which can be expected to have a significant positive effect not only on improving the accurate representation of the contribution of dynamical processes to SIC change in the model, but also on improving the simulation of sea ice thickness.