Abstract. In the last decade, various satellite missions have been monitoring the status of the cryosphere and its evolution. Besides sea-ice concentration data, available since the 1980s, sea-ice thickness retrievals are now ready to be used in global operational prediction and global reanalysis systems. Nevertheless, while univariate algorithms are commonly used to constrain sea-ice area or volume, multivariate approaches have not yet been employed due to the highly non-Gaussian distribution of sea-ice variables together with the low accuracy of thickness observations. This study extends a 3DVar system, called OceanVar, which is routinely employed in the production of global/regional operational/reanalysis products, to process sea-ice variables. The tangent/adjoint versions of an anamorphosis operator are used to locally transform the sea-ice anomalies into Gaussian control variables and back, minimizing in the latter space. The benefit achieved by such a transformation is described. Several sensitivity experiments are carried out using a suite of diverse datasets. The sole assimilation of the CryoSat-2 provides a good spatial representation of thickness distribution but still overestimates the total volume that requires the inclusion of Soil Moisture and Ocean Salinity (SMOS) mission data to converge towards the observation estimates. The intermittent availability of thickness data can lead to potential jumps in the evolution of the volume and requires a dedicated tuning. The use of the merged L4 product CS2SMOS shows the best skill score when validated against independent measurements during the melting season when satellite data are not available. This new sea-ice module is meant to simplify the future coupling with ocean variables.