Assimilation of ice compactness data in a strong coupling regime in the ocean - sea ice coupled model

. The Arctic Ocean plays an important role in the global climate system, where sea ice regulates the exchange of heat, moisture and momentum between the atmosphere and the ocean. A comprehensive analysis and forecast of the Arctic ocean system requires a detailed numerical ocean and sea ice coupled model supplemented by assimilation of observational data at appropriate time scales. 5 A new operative ocean – ice state forecast system was developed and implemented. It consists of the INMIO4.1 ocean general circulation model and the CICE5.1 sea ice dynamics and thermodynamics model with common spatial resolution of 0.25 ◦ . For the exchange of boundary conditions and service actions (data storage, time synchronization, etc.), the coupled model uses the Compact Modeling Framework (CMF3.0). Data assimilation is implemented in the form of the Data Assimilation Service (DAS) based on the Ensemble Optimal Interpolation (EnOI) method. This technique allows to simultaneously correct the 10 ocean (temperature, salinity, surface level) and ice (concentration) model ﬁelds in the DAS service, so they are coordinated not only through the exchange of boundary conditions, but already at the stage of data assimilation (i.e. strong coupling data assimilation). Experiments with the INMIO - CICE model show that the developed algorithm provides a signiﬁcant improvement in the accuracy of forecasting the state of the ice ﬁeld in the Arctic Ocean.

The practice of numerical weather prediction has shown that coupled models usually provide significantly more accurate forecasts than stand-alone ones. This has led to research on data assimilation (DA) in coupled models, which promises to further improve forecast quality (Skachko et al., 2019). An overview of current activities in coupled forecasting systems and coupled DA can be found in (Brassington et al., 2015). The World Meteorological Organization (WMO) Meeting on Coupled Data Assimilation (Penny and Hamill, 2017) defined the classification of weakly coupled and strongly coupled DA and their 5 variations.
The forecasting experience of leading scientific centers has shown that for an accurate mid-term forecast of the state of water and ice in the Arctic, the models of ocean and ice dynamics and thermodynamics need to be fully coupled. Thus, a DA technique capable of retaining sea ice in a coupled model in a dynamically consistent manner could help improving the accuracy of climate studies and weather predictions (Kimmritz et al., 2018). The goal of this work is to develop a coupled 10 ocean -sea ice model with coupled assimilation of available observational data and implement it with application for the Arctic region.
2 Coupled ocean-ice model with data assimilation

Coupled ocean-ice model setup
In this work, we use the global ocean -sea ice coupled model, which consists of the INMIO ocean general circulation model 15 (Ibrayev et al., 2012) and the CICE ice dynamics and thermodynamics model (Hunke et al., 2015) operating on massively parallel computers under control of the Compact Modeling Framework (CMF) (Kalmykov et al., 2018).
The ocean and sea-ice models use similar three-polar grids (Murray, 1996) with 0.25 • nominal resolution. The number of ocean z-grid horizons is 49 with vertical steps increasing from 6 m near the surface to 250 m in the deep part. The bottom topography data were interpolated from the (ETOPO5) array, excluding inland water bodies and small islands. 20 The INMIO-CICE model setup is the same as used in the ocean-ice component of (Fadeev et al., 2019). In particular, for the ocean model the lateral momentum exchange is modelled by the biharmonic operator with the coefficient equal to −1.5 · 10 11 m 4 /s on the equator and scaled proportionally to the grid cell area in power 3/2. The additional Smagorinsky biharmonic term in the form of (Griffies and Hallberg, 2000) is added to ensure numerical stability. The lateral heat and salt mixing is approximated by the Laplacian operator with coefficient 300m 2 /s on the equator and scaled proportionally to the 25 grid cell area in power 1/2. The momentum advection is approximated by the leap-frog scheme, while for the heat and salt advection the flux-corrected transport scheme (Zalesak, 1979) is used. The atmosphere-ocean surface fluxes of heat, mass and momentum are calculated by the (Launiainen and Vihma, 1990) atmospheric boundary layer bulk-formulae.
In the ice model configuration, the ice in each grid cell is split into five thickness categories with one additional category for snow. Thus, the array of ice concentrations (i.e. fractions of the grid cell area) for each category is used as the main variable 30 characterizing the state of the ice in the model. The elastic-viscous-plastic approximation is used to model the sea ice rheology, the ice transport is performed with the upwind scheme. To describe the evolution of temperature, a zero-layer thermodynamic model is used, in which ice is considered fresh and has zero heat capacity. A similar configuration of the CICE model is used in many key centers for operational forecasting of the ocean state, in particular in the TOPAZ system (Norway) consisting of a coupled ocean-sea ice model (HYCOM / CICE) with DA by the EnKF method (ensemble Kalman filter) (Sakov et al., 2012).
The time step for both INMIO and CICE components is 10 minutes, while their coupling interval is 20 minutes. The fields of upper grid layer temperature, salinity, velocity vector, freezing/melting potential, and surface tilt vector are sent from the ocean to the ice model. The fields (intergral over all categories) of ice concentration, horizontal ice-water stress vector, fluxes 5 of fresh water, salt, sensible heat and penetrating short-wave radiation are sent from the ice to the ocean model.
Before the experiments on the assimilation of observational data, the model spin-up was performed for the period of 2009.01.01 -2019.08.31, during which the atmospheric near-surface fields were defined by the ERA-Interim reanalysis (Dee et al., 2011) and the model solution was saved every 10 days. Further, these states, starting from the year 2011, were used as elements of the ensemble for the approximation of covariance matrices during the DA experiment (see Section 3). and December 1 2020. It can be seen from the figures that the model fields without assimilation are excessively smooth. No eddy dynamics is observed, which is typical for models with a coarse resolution. It is also worth noting that an overestimated amount of sea ice is produced by the coupled model, which is successfully corrected by assimilation. Also due to the assim-  ilation, the surface temperature field becomes more gradient, and the area of the ice cover decreases significantly, so that the model solution becomes more consistent with the OSTIA data.

EnOI Data Assimilation
The basic equations of the Ensemble Optimal Interpolation (EnOI) method are as follows (Evensen, 2003): where N is the ensemble size, and the matrix columns are equal to the model states minus their average over the ensemble.
Then the covariance matrix of model errors can be approximated on the basis of this ensemble in the following way:

Program implementation of EnOI
The assimilation procedure is encapsulated in the Data Assimilation Service (DAS) of the CMF framework. The service works in parallel on its separate CPU cores and can be used simultaneously for several model components (ocean, ice, atmosphere,  Table 1. Thus, assimilation of even one of the variables will correct the entire vector of model solution. Table 1. Contents and sizes of the vectors xa, xb, yobs, and the matrix A N b . T , S profiles data from the Argo project; absolute dynamic topography (ADT) along-track data from the Jason-3 satellite, AVISO project; sea ice concentration (SIC) data from the EUMETSAT OSI SAF project.

Element
Model state vectors (forecast xb and analysis xa) Observation data vector yobs In DAS, the seasonally varying ensemble (Oke et al., 2010) is implemented. For each call of the assimilation service, the model states for the same calendar month, but for previous years of model run, are used for the production of the ensemble.
Since only anomalies of the ensemble elements are used in the assimilation procedure, this approach allows to eliminate seasonal covariances, while preserving the mutual influence of synoptic circulation elements.

25
The following observational data will be used: So, the strongly coupled data assimilation approach was implemented in the ocean-ice model, when DA for different components is executed in the separate DAS service with coupled analysing fields from the two models. However, for smooth mutual adjustment, it was nevertheless required to implement an interface for the ice model, in which the correction factor field for ice concentration (K aice ) is introduced. The correction factor is equal to the model analysis obtained by the EnOI 15 algorithm (x a ) divided by the background (x b ), and is restricted to the range [1/α, α], where α is an empirical parameter. In our experiments, the optimal value of α is 4. At grid nodes where correction is not required or the model background has no ice the K aice coefficient is set to 1. Further, to correct the calculations in the CICE model, ice concentrations of various thickness categories (five ice categories and one snow category) are multiplied by the K aice field and the results are renormalized so that the integral concentration field is less than 1 at all nodes. 20 Thus, due to DA, there is a gradual "freezing" or "melting" of ice without direct correction of the ice thickness field.
The advantage over weakly coupled DA is that during the consequent correction of the ice concentration, other model fields (temperature, salinity, sea level) are also taken into account, thus they are coordinated not only through the mechanism of exchange of boundary conditions, but already at the stage of DA.
3 Experiments on simulation of water and sea ice circulation in the Arctic in 2020.
The proposed algorithm for data assimilation in a strongly coupled mode was tested in a series of numerical experiments with the global INMIO -CICE coupled ocean -ice model with 0.25 • resolution. The model was forced with atmospheric data of the 5 (GFS) forecasts accumulated for the period of 2020.01.01 -2020.12.31. Two experiments were carried out: 1. h01 -the control run without assimilation of observation data; 2. h02 -the experiment with strongly coupled assimilation of observed temperature, salinity, ocean surface level and sea ice concentration.
In this article, we will focus on the reproduction of the characteristics of sea ice only. Details of reproducing other charac-10 teristics of ocean circulation will be considered in separate articles.
3.1 Ice extent in the Arctic in 2020 Figure 2 and Table 2 compare the annaul variations of Northern Hemisphere sea ice extent in 2020 obtained in experiments h01 (no DAS) and h02 (DAS) with analysis data provided by the National Snow and Ice Data Center (NSIDC). It can be seen that due to the assimilation of observations, it was possible to significantly improve the accuracy of sea ice cover simulation.

15
The average of monthly errors decreased from 26 % to 6%. Figure 2 also shows that in the no-assimilation experiment h01, the resulting ice extent did not even fit into the "corridor" of interannual variability for 1979-2019. On the other hand, with assimilation in h02 we can see a clear correspondence of the simulation with the analysis data.  and after it (analysis). From these graphs we can conclude that assimilation gives the correct sign of the correction, and the difference between the RMS errors in the assimilating experiment and the control one is ∼ 0.5 • C. Furthermore, the absolute magnitude of the RMS error over all observation points for the SST is ∼ 1 • C, and ∼ 0.2 for the ice concentration. This is in agreement with data from the Copernicus Marine Environment Monitoring Service (CMEMS) (Melsom et al., 2019), in which 25 the RMS error of the daily forecast is ∼ 0.8 • C for SST , and 0.2 for the sea ice concentration. The methodology for calculating errors of sea ice concentration simulations is presented in (Desportes et al., 2017), (Melsom et al., 2019).

Conclusions
A strongly coupled data assimilation scheme was proposed and tested for the ocean -sea ice model INMIO-CICE with a wide range of observational data (Argo temperature and salinity, AVISO absolute dynamic topography of the sea level, OSISAF sea ice concentration). The main goal of the numerical experiments was to check the stable operation of the coupled iceocean prediction system, the correctness of the resulting model fields and the suitability of the entire system for operational