Toward a Data Assimilation System for Seamless Sea Ice Prediction Based on the AWI Climate Model

This paper describes and evaluates the assimilation component of a seamless sea ice prediction system, which is developed based on the fully coupled Alfred Wegener Institute, Helmholtz Center for Polar and Marine Research Climate Model (AWI‐CM, v1.1). Its ocean/ice component with unstructured‐mesh discretization and smoothly varying spatial resolution enables seamless sea ice prediction across a wide range of space and time scales. The model is complemented with the Parallel Data Assimilation Framework to assimilate observations in the ocean/ice component with an Ensemble Kalman Filter. The focus here is on the data assimilation of the prediction system. First, the performance of the system is tested in a perfect‐model setting with synthetic observations. The system exhibits no drift for multivariate assimilation, which is a prerequisite for the robustness of the system. Second, real observational data for sea ice concentration, thickness, drift, and sea surface temperature are assimilated. The analysis results are evaluated against independent in situ observations and reanalysis data. Further experiments that assimilate different combinations of variables are conducted to understand their individual impacts on the model state. In particular, assimilating sea ice drift improves the sea ice thickness estimate, and assimilating sea surface temperature is able to avert a circulation bias of the free‐running model in the Arctic Ocean at middepth. Finally, we present preliminary results obtained with an extended system where the atmosphere is constrained by nudging toward reanalysis data, revealing challenges that still need to be overcome to adapt the ocean/ice assimilation. We consider this system a prototype on the way toward strongly coupled data assimilation across all model components.


Introduction
Sea ice as one of the important components of the Earth System has come to prominence in recent decades owing to its dramatic retreat particularly in the Arctic (Kwok, 2018;Overland & Wang, 2010;Serreze et al., 2007). Sea ice prediction over a wide range of time scales is thus of great demand for both scientific and socioeconomic interests . A community collaborative effort was launched as the Sea Ice Prediction Network (SIPN) for both the Arctic and the Antarctic aiming to advance research on sea ice prediction. As part of SIPN since 2013, the Sea Ice Outlook (https://www.arcus.org/sipn/sea-ice-outlook) provides an open platform for researchers and stakeholders to share their summer ice prediction online (Blanchard-Wrigglesworth et al., 2017;Stroeve et al., 2014). Evaluation of the results from statistical models (e.g., linear regression models), sea ice-ocean coupled models, and fully coupled climate models shows that multimodel-mean sea ice prediction could generally be skillful for up to 5 months but remains challenging beyond this time scale (Wayand et al., 2019). However, real forecast skill still falls short of the inherent limits of sea ice predictability revealed by perfect-model type studies, (e.g., Day et al., 2016;Tietsche et al., 2014).

10.1029/2019MS001937
Key Points: • A data assimilation system is developed for seamless sea ice prediction on a climate model equipped with an unstructured ocean/ice component • The system performance is examined based on experiments assimilating synthetic (perfect model) and real observations • Multivariate assimilation of sea ice drift and sea surface temperature remarkably improves the ocean/ice model state in the polar regions Supporting Information: • Supporting Information S1 • Figure S1 • Figure S2 • Figure S3 • Figure S4 • Text S1 With numerous novel satellite sensors (e.g., Ice, Cloud,and land Elevation Satellite [ICESat-1/2], CryoSat-2, Soil Moisture Ocean Salinity [SMOS], and Soil Moisture Active Passive Mission [SMAP]) becoming available in recent years, sea ice thickness (SIT) retrievals and its further analysis over the entire Arctic are now possible. Ricker et al. (2017) merged the thickness data from the CryoSat-2 and SMOS into a weekly Arctic-wide CS2SMOS thickness data set using optimal interpolation (OI). Mu et al. (2018b) assimilated SIC together with SMOS and CryoSat-2 SIT, which are available only during the freezing season and thereby produced a thickness reanalysis (CMST) that also covers the melt season. Exploiting the brightness-temperature measurements from SMOS and SMAP, Paţilea et al. (2019) combined these two data sets with a better temporal and spatial coverage than SMOS thickness from Tian-Kunze et al. (2014). Impacts of the assimilation of these thickness data on sea ice simulations have been reported in several studies, (e.g., Allard et al., 2018;Lisaeter et al., 2007;Xie et al., 2016Xie et al., , 2018Yang et al., 2014. In the Southern Ocean, however, apart from the ASPeCt observations (Worby, 1999), the satellite SIT data are even more limited. ICESat SIT data only cover the period from 2004 to 2008 (Kurtz & Markus, 2012). Retrievals of SIT from CryoSat-2 are impeded by the insufficient snow knowledge in the Southern Ocean, and Environmental Satellite (ENVISAT) sea ice freeboard measurements are reported to have large biases (Schwegmann et al., 2016). Therefore, SIT estimates in the Antarctic are mostly provided by model output, such as Global Ice-Ocean Modeling and Assimilation System (Zhang & Rothrock, 2003), Southern Ocean State Estimate (Mazloff et al., 2010), and the model reconstruction by Massonnet et al. (2013).
Biases in the simulated seasonal SIT cycle have been found to cause wide sea ice spread in climate models . For sea ice prediction, it has been shown that better SIT initialization gives rise to significant improvement of sea ice forecasts from daily to seasonal scales in different systems. The importance of such long-term memory of SIT for sea ice prediction has been demonstrated in recent literature. Assimilation of thickness from SMOS (Yang et al., 2014) and CryoSat-2 (Mu, Yang, et al., 2018) improves thin-ice (0-1 m) and thick-ice forecasts, respectively, on daily scale in an ice-ocean coupled model (18 km). Using a high-resolution model (2 km) in the Beaufort Sea, Yaremchuk et al. (2019) also showed improvement of the 24-hr forecast skill. Utilizing thickness from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) to initialize a coupled ocean-atmosphere model, considerable changes in ice concentration forecasts have been demonstrated in the Coupled Forecast System model Version 2 (CFSv2, Saha et al., 2014) by Collow et al. (2015). Experiments in a perfect-model assimilation framework also suggest that winter thickness preconditions seasonal forecasts Zhang et al., 2018). Such reported advantages of winter SIT initialization can persist even longer. For example, Yang et al. (2019) used the ensemble initial conditions from the CMST system (Mu et al., 2018b) to conduct a seasonal forecast experiment from 25 May 2016, when direct satellite SIT was not available. However, promising SIC and SIT predictions in September were still found. The successful first trans-Arctic Passage of the Chinese icebreaker, supported by the same forecast system, also benefitted from the SIT assimilation in late April . Moreover, the atmospheric response to the improved SIT can also be considerable. For example, Blockley et al. (2018) found that SIT initialization in winter in Met-Office's Global Seasonal coupled prediction system further improves near-surface air temperature and pressure fields across the sea ice covered region.
Actual forecasts using standalone ocean/ice models not coupled to an atmospheric model component are limited by the boundary conditions. The climate models that evolve dynamically over time without prescribed boundary forcing are physically consistent across components. However, the systematic errors are easier to be amplified during long-term integration. Nevertheless, Sigmond et al. (2016) showed that dynamical forecast models outperform the simple statistical methods in their study. Zampieri et al. (2018Zampieri et al. ( , 2019 highlighted that the subseasonal ice prediction using coupled models can be skillful more than 1.5 months ahead in specific seasons. This leads to the question that whether a better initialization in all model components in a climate model particularly equipped with a high-resolution ocean configuration in polar regions could give rise to better sea ice prediction. Data assimilation approaches with relatively low computational cost such as nudging, OI, and 3D-var are widely used in climate models. The performance of these initialization methods varies across models due to different systematic biases or implementations (Wayand et al., 2019). The OI and 3D-var techniques need delicate treatment on the covariance matrix. The 4-D-variational method that is reported to be physically more consistent, however, calls for great efforts to implement in climate models. The EnKF dynamically estimates the covariance matrix and simplifies the implementation in complex models. For example, both the Norwegian Climate Prediction Model  and the Geophysical Fluid Dynamics Laboratory (Bushuk et al., 2017) apply EnKF data assimilation. Nevertheless, EnKF is still challenging when applied with highly nonlinear systems. In univariate assimilation, observations of specific variables are used to update only the same variables, which can lead to inconsistencies with other variables. In this case, adjustments of other model variables occur only indirectly through model physics and dynamics. By contrast, in multivariate data assimilation other variables in the same model compartment or even in other compartments ("strongly coupled assimilation") are updated based on cross-covariances (e.g., Massonnet et al., 2015;Sluka et al., 2016;Zhang et al., 2007). Multivariate assimilation is meant to reduce inconsistencies, but several challenges remain . These are related to, for example, dynamical inconsistences between the variables in the state vector during the coupled analysis; model biases that lead to systematic drift existing widely in climate models; and the imperfect approximation of the covariance matrix due to the insufficient sampling limited by the ensemble size. As a consequence, data assimilation in fully coupled climate models currently is conducted often in the perfect model framework in which model-simulated observations are assimilated, in other words, in the situation that the model has no biases (e.g., Kimmritz et al., 2018;Sluka et al., 2016;Zhang et al., 2007).
In this study, we describe and evaluate the assimilation component of a new seamless sea ice prediction system based on the AWI Climate Model (AWI-CM). The SIC, SIT, SID, and sea surface temperature (SST) observations are assimilated to initialize the model state. The seamless approach in this case refers not only to the prediction across a wide range of time scales but also to smoothly varying spatial resolution in the ocean-sea ice model. The latter is the unstructured-mesh Finite-Element Sea ice-Ocean Model (FESOM v1.4, https://fesom.de/; Sidorenko et al., 2011;Timmermann et al., 2009;Wang et al., 2008Wang et al., , 2014 developed at the Alfred Wegener Institute Helmholtz Center for Polar and Marine Research (AWI). It allows to conveniently use high resolution in the polar regions with the bulk of the global domain having coarser resolution. Consequently, the ocean and sea ice dynamics in high latitudes can be reasonably represented while the model computational cost remains relatively low. For data assimilation the system employs the Parallel Data Assimilation Framework (PDAF). Finally, we apply a simple nudging method to also constrain the atmospheric state in our coupled model. Schubert-Frisius et al. (2017) showed that spectral nudging of divergence and vorticity in atmosphere-only simulations can be used to downscale atmospheric reanalyses. In coupled models, it has been shown that the assimilation of atmospheric variables can improve ocean variables based on flow-depend covariances (Sluka et al., 2016).
The paper is structured as follows. The technical coupling of AWI-CM and PDAF is described in section 2. Section 3 provides the system design that aims at seamless sea ice prediction. In section 4, the system performance is evaluated, and the analysis increments are investigated. Discussions are given in section 5. Conclusions and perspectives are given in section 6.

AWI-CM
The AWI-CM (v1.1; Rackow et al., 2018;Sidorenko et al., 2015) consists of the ocean/ice model FESOM (v1.4) and the spectral atmosphere model ECHAM (v6.3.02p4; Stevens et al., 2013). FESOM is configured with variable resolution ranging from about 25 km in the northern North Atlantic and Arctic Ocean to nominal 1 • in most other parts of the global ocean. Using this model grid (CORE II; see global mesh on FESOM website https://fesom.de/models/meshessetups/ and Arctic mesh in Figure S1), it is possible to reasonably represent the Arctic and Southern Ocean sea ice in ice/ocean only simulations . We used the z coordinates option of FESOM with 47 unevenly spaced vertical levels. The model time step for FESOM is 900 s. A total of 288 cores are used for each FESOM run. For ECHAM we use the T63L47 configuration (up to 63 wave numbers and 47 vertical levels) with a time step of 400 s and 192 CPU cores for each run. The land surface model JSBACH (Stevens et al., 2013) and a hydrological discharge model (Hagemann & Dümenil, 1997) are included in ECHAM. FESOM and ECHAM exchange information every hour using the OASIS3-MCT coupler. The simulated mean climate state and climate variability exhibit good performance compared to established climate models that contributed to the fifth phase of the Coupled Model Intercomparison Project (CMIP5) as shown in Sidorenko et al. (2015) and Rackow et al. (2018).
The Finite-Element Sea Ice Model (FESIM)  is the sea ice component of FESOM. FESIM includes a single ice thickness category, zero-layer thermodynamics, and the elastic-viscous-plastic (EVP) solver for dynamics. Upon convergence, the EVP solver allows to reasonably represent small linear kinematic features in case of sufficient model resolution . For better numerical stability, the finite element flux corrected transport scheme is used to advect sea ice concentration and sea ice and snow thickness. The thermodynamic part follows Parkinson and Washington (1979) with a prognostic snow layer as in Owens and Lemke (1990). Snow to ice conversion takes place when the snow layer becomes thick enough to sink below the ocean surface. Both the sea ice and ocean are discretized on the same unstructured grid, so the fluxes between the two components can be exchanged directly. The atmosphere-ocean surface flux, including over sea ice, is calculated in ECHAM. The simulated spatial distribution of sea ice thickness is fairly realistic compared to other climate models as reported in Sidorenko et al. (2015).

PDAF
PDAF (http://pdaf.awi.de ;Nerger & Hiller, 2013) is an open source data assimilation framework designed to simplify the implementation of ensemble Kalman filters in numerical models in a computationally effective way. It currently provides filters like the Localized Singular Evolutive Interpolated Kalman (LSEIK; Pham et al., 1998;Pham, 2001) filter, the Local Ensemble Transform Kalman Filter (LETKF; Bishop et al., 2001), the Local Error Subspace Transform Kalman Filter (Nerger et al., 2012), and the original EnKF method (Evensen, 1994). The high efficiency and full parallelization features of PDAF allow for a large ensemble size (if necessary) for data assimilation in climate models. The online-coupled data assimilation circumvents the costly model restart during each assimilation cycle, which is crucial for climate models that usually need to request a huge amount of computation resources during job submission and write and read large amounts of restart data. In our current setup, the LESTKF in PDAF v1.13.1 is used. As described in Nerger et al. (2012), the LESTKF projects all the ensemble information onto the error subspace and then computes the ensemble transformation. This leads to higher computational efficiency of this filter compared to EnKF and LETKF and a better assimilation performance compared to LSEIK and EnKF.

Technical Implementation
The integration of AWI-CM and PDAF is illustrated in Figure 1. During the initialization, PDAF governs the adaption of the parallelization. Processors in the general communicator MPI_COMM_WORLD are grouped into different new defined communicators, COMM_filter, COMM_model, and COMM_couple. The COMM_filter provides processor communication for the filter. COMM_model carries out the communication among different climate model realizations. Communication between model realizations and the filter is performed in COMM_couple. This adaption on communicator endows AWI-CM with the capability to handle model simulations and data assimilation fully in parallel. The readers are referred to Nerger et al. (2019) for more details.
In the current setup, which has a focus on sea ice prediction, PDAF is implemented only for the sea ice and ocean model using LESTKF as shown in Figure 1. To additionally constrain the atmospheric model, the nudging technique for spectral variables can be applied in ECHAM. Strongly coupled assimilation with PDAF in both ocean and atmosphere components is planned to be explored for future versions (illustrated as the dashed box in Figure 1).
Twelve ensemble members and a localization radius of 126 km are used in the current setup. The chosen ensemble size is a compromise between computational costs and sampling of the covariances, noting that the relatively small localization radius avoids the use of often spurious large-distance covariances. The chosen localization radius is close to the zonally averaged first baroclinic Rossby radius of deformation around the latitude of 10 • (Chelton et al., 1998). Weights for the observations during analysis are calculated by the quasi-Gaussian weight function suggested by Gaspari and Cohn (1999). The AWI-CM ensemble members are run simultaneously on different cores. PDAF collects all the data directly from the memory through COMM_couple and conducts the analysis on the filter cores, then updates all the variables, writes diagnostic outputs, and distributes the analyzed variables back to the memory of the AWI-CM cores. The assimilation system restarts from the memory-stored analyzed fields, which speeds up the long-term integration.
Performing ensemble data assimilation with climate models is computationally demanding. Our configuration uses 5,760 processor cores. A typical wall time to perform an assimilation experiment of 1 year is around 2.25 hr. The exact cost depends strongly on the input/output (I/O), the computation platform, and how many variables are assimilated.

State Vector
The fields that are updated by data assimilation define a state vector. In the current assimilation system, the state vector consists of SIC, SIT, SID, and temperature and salinity in the mixed layer (Tmix and Smix). The mixed layer depth is determined at each assimilation step using the method of Levitus (1982) and is thus spatiotemporally variable. For a typical spatial distribution of the mixed layer depth, readers are referred to Figure 14 in Sidorenko et al. (2015). The omission of deeper-ocean temperature and salinity from the state vector reduces computational costs and is justified by earlier findings that the final performance of the ocean state analysis was not considerably improved in a similar assimilation system (Kimmritz et al., 2018). However, taking advantages of the cross-covariance between different variables, Tmix and Smix have been shown to benefit from ocean/ice assimilation in the surface (Kimmritz et al., 2018;Lisaeter et al., 2007).

Assimilation Data Sets
The real observation data used for assimilation are summarized in Table 1. The reprocessed OSI-450 SIC (Lavergne et al., 2019) from OSISAF is assimilated. The near-real-time SIC observation is used instead after the end of OSI-450 data on 1 January 2016. Although uncertainties are provided in the data set, a constant uncertainty of 0.15 is used in the system to account for the representation error (Janjić et al., 2017) and the satellite data error. Previous studies suggest that a constant uncertainty for SIC assimilation results in better summer SIC and SIT estimates . SIT estimates over the Arctic region based on CS2SMOS have been developed by AWI (Ricker et al., 2017). These weekly data are provided with a resolution of 25 km, which is close to the model resolution in the Arctic in our setup. It is thus well suited for data assimilation since the representation errors are expected to be small. To reduce abrupt changes that would occur if the data were assimilated weekly, the data are linearly interpolated to each day. In this case the uncertainties provided with the data set are used as the observation errors. Snow thickness assimilation is currently not performed due to its possible degradation on long-term sea ice extent estimates as reported in data with high quality defined by the status flag in the data set are used currently. This data set excludes data over the area with low ice concentration; thus, it has smaller coverage than other data sets. Therefore, smaller errors are introduced. Previous studies (e.g., Massonnet et al., 2014;Ungermann et al., 2017;Sumata et al., 2019) have already adopted this data set for data assimilation and parameter estimation.
SST from the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA; Stark et al., 2007) from the UKMO is used for the assimilation in the ocean. The OSTIA data are constructed by OI on a global 0.054 • grid with source data from multiple sensors and platforms (Table 1). The data record provides SST at UTC 00:00 and is thus free of diurnal variation, which is highly recommended by the data team for data assimilation in numerical weather prediction (NWP) systems. Several systems have already used this data set, such as TOPAZ4 (Sakov et al., 2012) and the ECMWF NWP system (Browne et al., 2019). Uncertainties provided in the data set are directly used as the observation error during assimilation.

Initial Construction of the Ensemble
The initial construction of the ensemble is performed by adding different perturbations to a single restart model state for each ensemble member. A reasonable choice of perturbations is required to ensure the stability of the model and data assimilation procedure. To obtain the perturbations, the SIC, SIT, SID, T, and S simulated in a free run from 1979 to 2013 are first collected. To ensure that the estimated probability space of the initial state is valid for the time of the year, from each year only the output over three consecutive days are used. For instance, in the situation that the system starts on 1 January, the model output on 31 December, 1 January, and 2 January are taken into account. Therefore, the matrix of state vectors consists of 105 columns. It is then decomposed into empirical orthogonal functions (EOFs) after subtracting its mean state. Second-order exact sampling (Pham, 2001) is applied to generate the initial perturbations. Specifically, the leading 11 EOFs are multiplied first with their corresponding singular values, second with a random matrix that preserves the mean and the covariances, and third with the square root of the ensemble size minus one. This ensemble represents the initial covariance matrix.

Preprocessing and Postprocessing
All the observations are preprocessed by removing unphysical values and data with low-quality flags. The Climate Data Operators (CDO) are used to interpolate from the observation grids to the unstructured model grid using distance-weighted average remapping. After the analysis, the absolute analysis increments of the variables in the state vector are limited to not exceed twice the ensemble spread (Sakov et al., 2012). Unphysical SIC or SST values generated by the filter are set to its physical bounds. For instance, when the analyzed SIC is lower than 0 or greater than 1, it is set to 0 or 1, respectively; the analyzed SST is limited to remain above the freezing point. These two procedures guarantee model stability and prevent excessive shocks during the initialization. Statistics from our experiments show that these limitations come into play only rarely. We find slight salinity drift when applying the real data assimilation due to systematic model errors. A correction of ∼10 −5 psu is added at the surface at every assimilation step to rectify the drift.

Experiment Design
To test the performance of the new system, we conduct the experiments listed in Table 2. To reduce computational costs, the different experiments are not run over the same time periods. The experiment CTRL is the reference historical run without assimilation. Exp_twin is the perfect model type experiment, in which synthetic SIC observations from CTRL are assimilated.
Exp_CTD_T assimilates real observations of SIC (C), SIT (T), and SID (D) in the sea ice model and SST (T) in the ocean model. In the experiment names we use the underscores to distinguish assimilation in different model components (sea ice model after the first underscore; ocean model after the second). For real SIC, SIT, and SST observations, they are assimilated daily; for SID, they are assimilated every other day. The initial conditions generated by Exp_CTD_T are further used as the restart files for the seamless sea ice prediction. In section 4.2, the analysis from this experiment is evaluated using different observations.
Three sensitivity experiments (Exp_C, Exp_C_T, and Exp_CTD_T_P) that assimilate different observations are also conducted (discussed in section 5). Particularly, Exp_CTD_T_P further constrains the atmosphere model state using reanalysis surface pressure data (SP, indicated with P in the experiment name) in ECHAM by nudging. To test the potential atmospheric effects on the assimilation results, the SP from ERA-Interim is nudged into the model. The original SP data are downscaled onto the ECHAM T63L47 grid and transformed into the spectral fields. The relaxation time scale is very short and carried out for all wave numbers and We evaluate the effect of the assimilation based on bias, RMSE, and ensemble spread. The ensemble mean states are used in the calculation of biases and RMSEs. Biases are calculated as X a − X re , where X a is the analysis and X re is the reference data. The RMSE is calculated as The spread is calculated as the standard deviation of the ensemble members. For all these metrics we compute area-weighted spatial averages over the area where the reference data are available.

Evaluation Data Sets
SIC observations retrieved from the Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) and the Advanced Microwave Scanning Radiometer 2 (AMSR2) with the ARTIST Sea Ice (ASI) algorithm (Spreen et al., 2008) are downloaded from the University of Bremen (https://seaice.uni-bremen. de/sea-ice-concentration/). We use the AMSR-E data with a resolution of 6.25 km, which cover the period from 1 June 2002 to 4 October 2011. AMSR2 data from summer 2012 to 2018 are used afterward due to the expiration of AMSR-E in October 2011. The AMSR-E/AMSR2 data can be considered as observations independent of the OSISAF SIC data, because the latter use a different retrieval algorithm (a combination of state-of-the-art algorithms and dynamic tie points) with multiple sensors on board different satellites.
ICESat SIT observations for the Arctic (Kwok et al., 2009) and Antarctic (Kurtz & Markus, 2012) in February and March 2008 are used for comparison. The uncertainty of ICESat SIT is quite large with 0.7 m in the Arctic and 0.23 m in the Antarctic. CMST (Mu et al., 2018a) and PIOMAS (v2.1; Zhang & Rothrock, 2003) data from 2011 to 2016 are downloaded and remapped to the EASE grid to facilitate the comparison with our assimilation results. These model-based data sets provide consecutive thickness estimates, but they may be also limited by the imperfect model dynamics; for example, CMST cannot reproduce the heavily deformed and ridged sea ice north of the Canadian Arctic Archipelago and Greenland (Mu et al., 2018a); PIOMAS tends to overestimate the thickness of thin ice and to underestimate the thickness of thick ice (Schweiger et al., 2011). Ice draft measured from upward looking sonar (ULS) moorings deployed in the Beaufort Gyre Exploration Project (BGEP; Krishfield et al., 2014) is converted into ice thickness by multiplying with a factor of 1.1 as suggested in Rothrock et al. (2008). The IceBridge SIT data (Kurtz et al., 2013) from both the calibrated product (Version 2.0) and the Quicklook product are collected. The SIT at the nearest model grid points are selected to compare with the IceBridge data.
Quality-controlled subsurface ocean temperature and salinity profiles from the objective-analysis EN4 data set of the UKMO (Good et al., 2013) are used for assessing the hydrography. For evaluating the atmospheric states, 10 m wind and 2 m air temperature reanalysis data from ERA-Interim are used.

Results
In section 4.1 we evaluate the twin experiment in which the output from the control run is considered as the truth and is assimilated back into the system. The real data assimilation is then described in Figure 2. Bias (a,b,c,d) and RMSE (e,f,g,h) of sea ice concentration (SIC) and sea ice thickness (SIT) in the Arctic (left column) and the Antarctic (right column) for CTRL (blue) and Exp_twin (orange). The unit for SIT is meter. SIC ranges from 0 to 1. X-axis represents years from the start. section 4.2, in which the analysis from Exp_CTD_T is further compared with independent observations and reanalysis data.

Synthetic-Data Assimilation
We conduct this twin experiment for two main reasons. The first reason is to assess the fundamental ability of the system to assimilate data without undesired side effects. The second reason is to diagnose whether multivariate assimilation leads to artificial drift during long-term integrations in the absence of model biases.
Synthetic SIC observations are constructed from the historical run (CTRL in section 3.5) of AWI-CM produced for the sixth phase of the Coupled Model Intercomparison Project (CMIP6), using the period from 1997 to 2004. Like for the real SIC data, observation errors of 0.15 are used for these synthetic data. The synthetic SIC observations are assimilated once per day, which is different from the monthly assimilation in Kimmritz et al. (2018). Detrimental system drift, if present, is expected to be amplified with higher assimilation frequency. The multivariate assimilation also updates SIT that is typically positively correlated with SIC (Lisaeter et al., 2007;Tietsche et al., 2013). Therefore, stronger drift would be expected for SIT than for other variables that have weaker covariance with SIC, such as sea surface salinity (SSS), as reported in Lisaeter et al. (2007) and also found in our experiments.
We show the biases and RMSEs calculated with respect to the truth for the CTRL and Exp_twin experiments in Figure 2. All the biases and RMSEs in the Arctic and Antarctic are reduced due to the SIC assimilation in Exp_twin. SIC biases and RMSEs in Exp_twin are less than 0.03, which suggests the SIC observations are assimilated properly. SIC assimilation does not pull the SIT state exactly back to the truth; however, bias and RMSE reductions of nearly 60% are found for SIT in Exp_twin. No apparent drifts for SIT biases and RMSEs are observed (Figures 2c, 2d, 2g, and 2h). This indicates that the ensemble sampling and localization radius do not introduce significant artificial degradation of the state. However, the obtained seasonality in temporal evolution of SIT bias indicates seasonal weakening in the sampled and represented correlations between SIC and SIT. Overall, this experiment shows that if the model is unbiased and assimilated data are perfect, the assimilation system works properly in the current configuration.

Real-Data Assimilation 4.2.1. Sea Ice 4.2.1.1. Sea ice Concentration
We compare the SIC analysis with the independent AMSR-E/AMSR2 data in Figure 3a. The RMSE of SIC is strongly reduced both in the Arctic and Antarctic, as expected. The RMSE of SIC for Exp_CTD_T is 0.17 and 0.19 in the Arctic and Antarctic, respectively. This is quite close to the OSISAF data (0.16 and 0.18). An annual cycle is found in the RMSE evolution. In the melt season, RMSEs are larger than in the freezing season in both hemispheres. The Arctic RMSE is smaller than that of the Antarctic. RMSEs increase starting from 2016, likely due to the change of the OSISAF SIC product used for assimilation. This might be related to the quality of the data since a slight increase of the OSISAF RMSE is also observed at the beginning of 2016 (Figure 3a). Furthermore, using the same observation error for assimilating the new data may not be optimal.
The maximum ensemble spread of SIC each year does not change much before 2016 in both the Arctic and the Antarctic (Figure 3b). Annual variability is also found in the evolution of the ensemble spread, which is larger in winter and smaller in summer. Smaller spread in the freezing season is expected because of the saturation of SIC (100%). A smaller ensemble spread of SIC corresponds to a larger RMSE of SIC. Such correspondence suggests the seasonality of the model background (forecast) error. The larger ensemble spread and RMSE in the Antarctic compared to the Arctic Ocean indicates that the ocean drives more sea ice variations in the Southern Ocean than in the semi-closed Arctic Ocean.

Sea Ice Thickness
For the sea ice thickness evaluation, there are no independent long-term and Arctic-wide sea ice thickness records. Thus, we first consider the CS2SMOS data that were also assimilated. Exp_CTD_T has an averaged SIT RMSE of 0.19 m. This is a clear reduction compared with a value of 0.63 in the CTRL (Figure 4a) experiment. Changing the SIC product for assimilation in 2016 leads to a larger SIT RMSE as also found in SIC RMSE. The ensemble spread of SIT, however, is reduced after the assimilation of SIT starting from October 2010 (Figure 4b). Changing the SIC product introduces no obvious systematic effects on the SIT spread, although values in 2018 are unusual.
The SIT distributions in February and March are further compared with ICESat and CMST data ( Figure 5). In 2008, there is no SIT assimilation in both hemispheres, so differences between CTRL and Exp_CTD_T are solely due to the assimilation of other variables through nonzero cross-covariances and model dynamics. In  the Arctic, the SIT derived from ICESat is reported to be thinner than ULS observations in the Beaufort Sea (Kwok & Cunningham, 2008). This indicates that the larger thickness in the model simulations compared to ICESat might be reasonable (Figures 5a and 5d). As a result, we cannot conclude whether the assimilation improves the thickness estimate in the Arctic Ocean from this comparison. However, the comparison with the CMST data reveals that the mean March SIT over 2011-2016 in Exp_CTD_T is notably improved over perennial ice areas in the Arctic, which is underestimated strongly (more than 1.0 m) in the CTRL experiment (Figures 5c and 5f). Assimilation of SIC and SID in the sea ice model gives rise to SIT increments in the Antarctic in February and March (Figures 5b and 5e). The assimilation brings the SIT closer to the ICE-Sat data in the Southern Ocean even without any SIT assimilation. The ice thickness in the ICESat data has been shown to be around 0.23 m thicker than according to the ASPeCt observations during that time (Kurtz & Markus, 2012). Even so, the reduced negative thickness bias in Exp_CTD_T versus CTRL relative to ICE-Sat suggests an improvement of the model state. Considering the period 2011-2016 for which the CMST data are available and during which CS2SMOS have been assimilated, the mean March SIT in Exp_CTD_T shows notable improvement over perennial ice areas in the Arctic. The thickness is underestimated substantially (more than 1.0 m) north of Greenland in the CTRL experiment.
We now consider September CMST and PIOMAS data to evaluate the assimilation in late summer in the Arctic Ocean. Bias reductions are apparent with respect to both data sets ( Figure 6). Differences between the analyzed SIT and CMST are rather small. The reason is that the two systems assimilated thickness data from the same sources, although the approaches are quite different (Mu et al., 2018b). Compared to PIOMAS, the SIT analysis from Exp_CTD_T is thicker particularly north of the Canadian Arctic Archipelago (more than 1.0 m). In the Arctic interior, the SIT analysis is generally 0.3 m thicker than PIOMAS. Note that PIOMAS appears to underestimate the SIT over thick ice areas as documented in earlier studies (e.g., Mu et al., 2018b;Schweiger et al., 2011).
The seasonal cycle for the SIT analysis is compared with three ULS moorings in the Beaufort Gyre (Figure 7). The CTRL experiment overestimates SIT at all three locations as shown by both the plots of seasonal variations and the scatter plots of SIT. The assimilation significantly improves the representation of the SIT. The mean SIT analysis at the three locations is 1.0, 1.3, and 1.2 m, respectively, which is close to the ULS observations with 1.1, 1.3, and 1.3 m. We find underestimation in spring and overestimation in winter in the SIT analysis compared to the ULS observations. This can be explained by the differences between the assimilated CS2SMOS SIT data and the ULS data (Figures 7a-7c).
The IceBridge data from 2009 to 2018 are used to assess the simulated SIT further. Scatter plots (Figure 8) illustrate that the SIT analysis fits better to the observations than the CTRL run. The RMSEs with respect to IceBridge data for CTRL and Exp_CTD_T are 0.95 and 0.80 m, respectively. The coefficient of determination (R 2 ) is 0.21 for CTRL and 0.47 for Exp_CTD_T. Thin SIT in the 1.0-2.0 m bin is overestimated in the CTRL run. This bias is reduced with data assimilation. However, the SIT analysis still underestimates the observed thickness, especially when IceBridge thickness is over 4.0 m.  (d-f) The same as (a)-(c) but for the Antarctic. For the Arctic, the mean March SID is shown, while the mean September SID is shown for the Antarctic. Note that only the data with high quality are plotted for OSISAF.

Sea Ice Drift
The Arctic SID simulated in CTRL over 2008-2018 exhibits a pattern representing a climate state (positive Arctic Oscillation) opposite to that in the OSISAF observations (negative Arctic Oscillation; cf. Figures 9a  and 9c). On the one hand, the more realistically represented SID in Exp_CTD_T (Figure 9b) can be due to the assimilation of OSISAF drift data directly. However, it can also be an indirect effect caused by the more realistic SST pattern, which considerably improves the Arctic sea level pressure distribution and atmospheric circulation, which in turn leads to a better SID pattern. The Antarctic SID pattern is simulated quite well already in the CTRL run (Figure 9d). Changes in the SID after applying the assimilation are found in particular along the Antarctic coastline, where the SID analysis has a lower ice speed (Figure 9e). The effects of the reduction in sea ice speed are further discussed in section 5.1.

Ocean
We calculate the RMSE of the ocean temperature at 5 m depth relative to the EN4 analysis data over the assimilation period. The global mean temperature RMSE of the CTRL experiment is 1.78 • C. This error is reduced by ∼70% in the Exp_CTD_T experiment. As shown by the spatial SST differences in Figures 10a  and 10b, the systematic errors in AWI-CM reported by Sidorenko et al. (2015) are considerably reduced, as expected. It is worth noting that over the areas with westward-intensified currents (e.g., Gulf Stream, Kuroshio Current, and Brazil Current), our SST analysis still shows some weak biases (Figure 10b). This suggests that only applying SST assimilation is not sufficient to constrain the model temperature in these current-separation areas. Sea surface height assimilation is expected to correct such biases in future versions.
Obvious improvements of temperature are also found at 50 m depth. Most biases are eliminated in the analysis in comparison to the CTRL run (Figures 10c and 10d). Besides the western boundary current areas, some biases still exist over the east tropical Pacific Ocean, the tropical Indian Ocean, and the upwelling area near western Africa in the analysis (Figure 10d). However, they are much smaller in magnitude compared to the CTRL run. Some warm biases (∼2.0 • C) are induced over the tropical Indian Ocean.
The temperature patterns at the depth of 150 m (Figures 10e and 10f) and 300 m (Figures 10g and 10h) suggest that there are still bias reductions, but much smaller than those in the upper ocean. The temperature analysis in the deeper Indian Ocean is deteriorated by the assimilation. We suspect that the positive SST correction enhances the stratification in the vertical column over this area, which impedes the oceanic heat release. The vertical profile of temperature in Figure 11a shows that the improvements can penetrate to a depth of about 300 m. The deep ocean (400-1,000 m) is also improved, likely through dynamic processes, but the improvement is minor (Figure 11a).
The salinity bias relative to the EN4 data exhibits limited improvement in Exp_CTD_T compared to that in CTRL ( Figures S2 and 11b). The maximum improvement in the vertical salinity profile is about 0.028 psu, which is about 14% of the total error. We used the covariance between SST and salinity to also update the salinity in the mixed layer (Smix). The relationship between SST and Smix, however, is not always consistent even in perfect-model experiments, as documented in Zhang et al. (2007). This further gives rise to varied salinity changes in different ocean basins ( Figure S2). For example, the assimilation introduces negative biases northeast of Australia ( Figures S2b and S2d), while surprising bias reductions are found over the tropical Atlantic Ocean (Figures S2a and S2b). Systematic model bias, insufficient ensemble size, and nonlinear dynamics can lead to inconsistencies when multivariate assimilation is applied. The assimilation of salinity profiles holds promise to prevent spurious salinity states, which can be crucial for sea ice prediction in high latitudes where salinity plays an important role for the ocean circulation.

Atmosphere
Through ocean-atmosphere interactions, the atmospheric state is influenced when SST observations are assimilated. Not surprisingly, the bias of 2 m temperature over the ocean, which generally reflects the SST bias shown in Figure 10a, is systematically reduced (cf. Figures 12a and 12b). The biases over the continents Figure 11. Lateral mean profiles of potential temperature (a) and salinity (b) averaged over the ice-free area. The biases of the vertical profiles relative to the EN4 data are also shown in the inset plots. largely remain, although some improvements can also be seen over land, for example, in South America. These results are expected as 2 m temperature is strongly influenced by the underlying surface. Notably, the temperature over the Arctic Ocean is also slightly (∼0.5 • C) improved (Figure 12b) directly because of the SST assimilation in the melt season, and possibly also indirectly because of the improved ocean and atmosphere states in lower latitudes.
Winds are influenced by SST changes less directly than near-surface temperatures. Nevertheless, errors in the zonal 10 m wind component over the tropical oceans in Exp_CTD_T are considerably reduced (Figures 12c and 12d). The biases of the meridional 10 m wind component are reduced even more over the ocean, in particular in the region of the Intertropical Convergence Zone (Figures 12e and 12f). These give rise to a positive feedback over the upwelling area. For instance, the SST bias in the tropical ocean east of Africa is possibly caused by the insufficient wind stress due to the coarse atmosphere model. The northward meridional wind over this region is weaker than in the reanalysis data, which implies a weaker upwelling of cold water. The SST assimilation brings about more realistic temperature and sea level pressure distributions that induce stronger wind, which is in closer agreement with the reanalysis (Figure 12f). The stronger wind lifts more cold water, cools down the surface, and further reduces the SST bias (Figures 10a and 10b).

Analysis Increments in the Exp_CTD_T
Mean analysis increments reveal systematic model errors and thereby provide hints for further model development. For sea ice, large mean SIC (>0.005) and SIT (>0.01 m) analysis increments mostly occur in the Marginal Ice Zone in the freezing season for both hemispheres, while they are found over the entire ice-covered area in the melt season (Figure not shown). The seasonal variations of the spatially averaged analysis increments are more instructive. The SIC and SIT increments in the Antarctic are quite consistent as both are positive from November to March and negative from April to October (Figure 13). In the Arctic, the seasonal cycle has a phase opposite to that in the Southern Ocean. The conspicuous negative SIT increment in October (orange dotted line in Figure 13) can be attributed to the sudden onset of the SIT assimilation when CS2SMOS data become available in October. Nevertheless, there is largely a positive relationship between SIC and SIT increments in the Arctic Ocean, too.
The increments show that the model is seasonally biased. The model generally overestimates sea ice melt in the melting season and overestimates sea ice freezing in the freezing season ( Figure 13) and thus overestimates the seasonal cycle. Even with SIT assimilation in the Arctic, the SIT is still overestimated in the freezing season. Note that there is no SIT assimilation in the Antarctic but the increments there generally deliver the same message as in the Arctic. The results indicate that some physical processes in the sea ice Figure 12. Biases of the mean 2 m temperature (a and b), zonal (c and d), and meridional 10 m (e and f) wind components relative to ERA-Interim data averaged over the assimilation period. The left column is for CTRL, and the right is for Exp_CTD_T. model might need improvement through parameter estimation or more fundamental changes. For example, FESOM uses a zero-layer approximation where the sea ice heat capacity is neglected. With heat capacity, sea ice freezing can release more heat, while sea ice melting needs more heat (Semtner, 1976). With the same atmosphere and ocean forcing, this will cause sea ice to grow slower in winter and melt slower in summer, which can compensate the seasonal bias. A new sea ice model version accounting for heat capacity is under development and will be used in future versions of the system. The seasonal cycle of analysis increments can also be used for bias corrections prior to the analysis step in the assimilation system.
The mean SID increments in winter for both hemispheres are shown in Figure S3. The SID increments in March reveal that the assimilation causes slower drift along the East Greenland Current (EGC), faster drift in the Beaufort Gyre, and more onshore drift in the Siberian Seas ( Figure S3a). The most pronounced effect in the Antarctic is the deceleration of the ice drift along the Antarctic coast ( Figure S3b) that is also found in the mean state (Figure 9). The possible reasons, however, are difficult to disentangle in a coupled climate model where different model components may play a role for biases and thus the analysis increments. For example, in the Fram Strait and downstream of the EGC, where sea ice rheology plays less important role, too strong northerly winds simulated by the atmosphere model can lead to the bias there. In the Southern Ocean, the slowdown of SID along the Antarctic coast in the analysis may be due to the incapability of the sea ice model to simulate the ridged thick ice in those regions ( Figure S3b). In addition, the coarse atmosphere model inevitably introduces errors due to the interpolation from the atmospheric grid to the ocean grid at the coast. To solve these model deficiencies, one possible solution could be to optimize the ice-ocean exchange coefficients and the ice-atmosphere exchange coefficients used in the interfacial stress calculation (Massonnet et al., 2014).

Assimilating SID Improves SIT Estimate
The twin experiment Exp_twin that reduces the model bias by assimilating the synthetic observations exhibits virtually no system drift ( Figure 2). However, drift is evident when it comes to the real data assimilation. For example, when only assimilating SIC but also updating SIT by exploiting their covariance matrix in Exp_C_T, the Antarctic SIT along the coast reaches more than 8.0 m (Figure 14a) starting from initial values of ∼2.0 m. Ridged thick ice along the Antarctic coast can occur in principle but has not been observed extending over such a large area. The same situation is also found in Exp_C (Figure not shown) indicating that this behavior is linked to systematic errors introduced by the SIC assimilation, not the SST assimilation.
When SID is assimilated into the system in Exp_CTD_T, more realistic SIT around the Antarctic is achieved (Figure 14b). Note that there is no SIT assimilation in the Antarctic in Exp_CTD_T because no observational thickness data are available. The SID generally has short-term memory because of large internal stress in sea ice. The atmosphere, ocean, and internal stresses quickly undo any changes to the ice drift pattern immediately after the analysis. However, the result shows that the SID assimilation does nevertheless help to obtain  a better SIT analysis. A previous study on the Arctic Ocean also found that not only the SID itself but also the SIT and the SIC over pack ice area can be improved with SID assimilation ).
The assimilation of SID reduces the offshore drift that could give rise to sea ice accumulation (Figure 9), thus resulting in thinner sea ice along the coast. We find that the covariance between SIT and SID has magnitudes similar to those between SIT and SIC. However, when taking the innovation and observation errors into account, the weight given to the SID assimilation is about three times larger than that given to the SIC assimilation. That is, the SID assimilation plays a more important role for updating SIT during the analysis. The tendency for SIT to become too large induced by the systematic SIC errors is compensated by the improvement resulting from assimilating SID. Therefore, the SIT estimate can be improved through assimilating SID even when SIT observations are not available. It also implies that we should use as many multivariate data as possible for state estimates when observations of some of the variables are not available. Apparently, these data should bring about realistic cross-covariances and sufficiently small biases. This also resolves the apparent contradiction that SID has a significant impact despite the negligible memory of SID itself.

SST Assimilation Corrects the Arctic Circumpolar Boundary Current
The Arctic Circumpolar Boundary Current (ACBC) is a counterclockwise circulation of water masses (Atlantic Water [AW]) following the topography as shown by observations and high-resolution models (Aksenov et al., 2011;Wang et al., 2018). However, a reversed current has often been simulated in coarse sea ice-ocean models (Holloway et al., 2007). Based on idealized, singlelayer and barotropic model experiments, Yang (2005) suggested that the potential vorticity (PV) budget determines the direction of the boundary current. Analyses based on coordinated simulations (Karcher et al., 2007) found that the net flux of PV from the Barents Sea controls the circulation mostly in the Eurasian Basin, while the vertical PV flux (wind driven) is more important for the Beaufort Sea. In the theoretical and model study by Spall (2013), the sensitivity of this circulation to model parameters (e.g., vertical mixing coefficient) is explored, and the balance between the barotropic effect (wind driven sea surface height gradient) and the baroclinic effect (sloping isohalines) is proposed to determine the circulation. Despite improved theoretical understanding, it is still a challenge to simulate a realistic AW layer in coupled climate models with generally coarse resolution (Shu et al., 2019). In our climate model CTRL run, although a relatively high resolution (∼25 km in the Arctic) is used, a reversed current (clockwise) is simulated in the Canadian Basin as shown by the plot of topostrophy (Figure 15a). Topostrophy was introduced by Holloway (2008) to transform velocity vectors to scalars. A positive topostrophy corresponds to flow that has deep water to the left on the Northern Hemisphere. In the Arctic Ocean, it then represents a counterclockwise flow along the boundary, and vice versa. Since the stand-alone ocean model in AWI-CM forced by atmospheric reanalysis data gives correct (counterclockwise) ACBC in the Arctic Ocean (Wang et al., 2018), the problem in the coupled climate model can be attributed primarily to the atmosphere component.
The atmosphere model ECHAM has a resolution of T63 (∼2 • ) in our simulations. As in most climate models, the sea level pressure in the historical run of the AWI-CM has higher pressure (∼3 hpa) than ERA reanalysis data over the northeast Barents Sea and Kara Sea. Over the Beaufort Sea, the sea level pressure is lower by ∼5 hpa (figure not shown). These biases thus reduce the pressure gradient and lead to an extension of the Beaufort Height further toward the Eurasian Basin. The bias in surface winds then influences the AW layer circulation through impacting the barotropic and baroclinic balance as discussed by Spall (2013).
The assimilation of SST in climate models can provide a better surface boundary condition for both the ocean and atmosphere. The improvement of atmosphere states over the ice-free area can be further propagated to the whole Arctic dynamically. Such improvement is only found in Exp_CTD_T and Exp_C_T, but not in Exp_C (figure not shown), indicating that the improvement is due to the SST assimilation. The more realistic atmospheric circulation appears to correct the direction of the ACBC in the Arctic Ocean as shown by the topostrophy (Figure 15b). No considerable temperature and salinity improvements are found deeper than 150 m over the AW branch along the Norwegian shelf in Exp_CTD_T (Figures 10 and S2). This further confirms that the circulation change in the Arctic AW layer is associated with the change of the atmospheric circulation. As the atmospheric circulation in Exp_C is not improved, the boundary current is still in the wrong direction over the whole integration period (Figure not shown).
We note that the completely corrected circulation direction in Exp_CTD_T occurs starting from 2015. So, it takes 8 years for the assimilation to reverse the clockwise circulation to the realistic counterclockwise circulation. This suggests that assimilation over a certain time period is required to reconstruct the baroclinic structure of the Arctic AW layer.

Toward a Better Assimilation System Design
The performance of EnKF depends considerably on the ensemble size, the localization radius, and the spread. The ensemble size we used here is 12, which may be not large enough to compute a realistic covariance matrix. However, we did test a higher ensemble size of 24. The larger ensemble does improve the analysis, but not to the extent that it would justify the increase of computational cost by a factor of 2. The localization radius we used is constant in space. It could also be set to a function of the local first baroclinic Rossby radius to optimally exploit the information contained in ensemble-estimated covariances, which has not yet been explored.
In the setups discussed above, the ensemble spread is maintained at a sufficiently high level because the atmospheric trajectories become quite different from each other after a short period of integration; the ocean and ice of each ensemble member essentially experience their own weather-type atmospheric variability. For seamless sea ice prediction that is meaningful also at short lead times, it would be important to constrain the initial atmosphere states, too. Furthermore, it is expected that better constraining the atmospheric state should help to initialize the ocean and sea ice. To explore this, we conducted a 1-year experiment Exp_CTD_T_P where sea surface pressure data from the ERA Interim reanalysis are nudged into the system. The atmosphere is constrained as expected. The improvement in the atmosphere is observed not only at the surface but also at higher levels of the atmosphere. The sea ice and ocean analyses, however, are deteriorated with a SIC RMSE increase of 0.08 in the Arctic and 0.09 in the Antarctic with respect to Exp_CTD_T. The likely reason is that the ensemble spread is decreased due to the similarity in the atmosphere states. In consequence, the model background receives more weight compared to the observations. Additional ensemble inflation and/or further tuning are thus required when we implement atmospheric data assimilation in the future.
To reduce systematic biases in real data assimilation, there are several approaches, such as forecast correction (Barth et al., 2015), bias propagation (e.g., Nerger & Gregg, 2008), and parameter estimation (e.g., Annan et al., 2005). We also did an experiment that has only SIC and SIT in the state vector with only SIC assimilation. Because of the systematic SIC bias, the SIT exhibits a drift in the 10-year integration. We fit the analysis increment of SIC to SIT and get an exponential function for each month. This function is then used to deduce the background bias in each assimilation step. The SIT drift caused by the multivariate assimilation due to the systematic bias from assimilating SIC is removed, and the SIT field is better estimated. This promising preliminary result suggests that such approaches merit further research in the context of multivariate data assimilation.
We aim to assimilate more observations because systematic errors can be greatly reduced if observational data are available. Sea surface height data from satellite observations can impact ocean currents significantly, while temperature and salinity (T/S) profiles can provide long-term memory for seasonal forecasts.

10.1029/2019MS001937
These data are not assimilated currently. Moreover, the inconsistency between different variables used in the assimilation and between different model components during the assimilation is still challenging. Further related research is required to explore, for example, whether the simulated covariance between SIC and salinity in polar regions is realistic and to what extent the non-Gaussian distribution of sea ice can affect other variables during the assimilation.

Conclusions and Perspectives
High-quality initialization is crucial for climate model-based prediction. Therefore, data assimilation is at the center of any operational forecast system. The ensemble-based data assimilation method is easier to be applied in complicated climate models compared to the variational approach. In this study, the prototype of a new seamless sea ice prediction system based on AWI-CM1.1 is described, and the performance of the assimilation is evaluated. A twin experiment assimilating synthetic observations shows robustness of the assimilation system. In the real data assimilation, the system assimilates SIC, SIT, SID, and SST. The evaluation against in situ and reanalysis data indicates that the system performs well. Using experiments assimilating different combinations of observations, we found that the assimilation of SID reduces the systematic bias induced by assimilating SIC and improves SIT estimates substantially. We further show that assimilating SST in ice-free regions improves the overall ocean-ice-atmosphere system and even corrects the direction of the Atlantic Water circulation at depth in the Arctic Ocean.
Moreover, AWI-CM has been upgraded to the second version (Sidorenko et al., 2019) recently with a model speedup of approximately five times compared to the first version in terms of computational efficiency of the ocean/ice component. Using the new version will enable us to use a larger ensemble size, higher resolution, and/or more observations such as T/S profiles and sea surface height. In particular, using higher model resolution will not only better simulate sea ice leads  but also improve the representation of landfast ice, both of which are important for high-latitude climate. We also plan to further extend the data assimilation system with PDAF to the atmosphere model and thereby to conduct consistent (strongly coupled) data assimilation linking different model components. The primary goal of the continuing development is to use multivariate and fully coupled assimilation in AWI-CM to provide skillful and efficient sea ice prediction. The assimilation system presented here paves the way to tackle actual forecasts beyond the assimilation phase in the near future.