Arctic Sea Ice Thickness Estimation Based on CryoSat-2 Radar Altimeter and Sentinel-1 Dual-Polarized SAR

We present a method to combine CryoSat-2 (CS-2) radar altimeter and Sentinel-1 synthetic aperture radar (SAR) data to obtain sea ice thickness (SIT) estimates for the Barents and Kara Seas. Our approach yields larger spatial coverage and better accuracy compared to estimates based on either CS-2 or SAR alone. The SIT estimation method developed here is based on interpolation and extrapolation of CS-2 sea ice thickness (SIT) utilizing SAR segmentation and segmentwise SAR texture features. The SIT results are compared to SIT data derived from the AARI ice charts, to ORAS5. PIOMAS and TOPAZ4 5 ocean-sea ice data assimilation system reanalyses, and to daily MODIS based ice thickness charts. Our results are directly applicable to the future CRISTAL mission and Copenicus programme SAR missions.

• C to represent dry snow cold conditions, and T am >0 • C to moist/wet snow conditions where the snow pack, if exists, prevents or significantly attenuates radar returns from sea ice. Between these T am limits the state of the ice pack is ambiguous as it depends on T a history and time of the CS-2 and S-1 data acquisition, e.g. early morning vs. late afternoon. The T am data will not be used classify the CS-2 and S-1 into different weather condition classes, but to support analysis on the accuracy of SIT 95 from the combined CS-2 and S-1 data. It is noted that the T am limits above are only approximations.
At the Vize station, representing the northern part of the Kara Sea, T am during Janury-April 2016 is mostly below -5 • C, there is only few instances of one to three day long periods of T am above -5 • C, and −T am is all the time below 0 • C. T am at the Popova station, which is located at the central Kara Sea, has similar behaviour, but from around 21 Apr onwards T am is over -5 • C. At the Kongsøya station located on the Barents Sea, east of Svalbard, T am is much warmer and is below -5 • C 100 only during half of the days in January-April, but still T am is over 0 • C only in the beginning of Jan. At the Varandey station in the Pechora Sea T am shows very cold weather in Jan, short periods of T am >-5 • C in Feb and Mar, and then rather warm weather after first week of Apr. We estimate that the period Jan-Apr 2016 represent typically cold dry snow conditions in our study area.
During October 2016 to Apr 2017 T am from the Vize station is over 0 • C during the first week of October and mostly over 105 -5 • C up to end of third week of November. After that there is only few cases of 1-3 day long periods of T am above -5 • C. At the Popova station T am is over 0 • C during first two weeks of October, and mostly below -5 • C after beginning of November.
The Kongsøya T am is all the time over -5 • C up to 24 November, and afterwards there are variable periods of T am above and below -5 • C with an average decreasing T am trend. T am is all the time below -5 • C in Apr. At the Varandey station T am is above 0 • C up to end of October. After beginning of November T am is mostly below -5 • C. There is a period of warm weather 110 in 15-26 Mar with T am sometimes being slightly over 0 • C. We conclude that CS-1 and S-1 data from the beginning of October to roughly mid-November is typically affected by wet or moist snow conditions, depending on the diurnal acquisition time and location in our study area (more probable in the Barents and Pechora Seas). Data acquired from mid-November 2016 to Apr 2017 mostly represent cold winter conditions.
In October -December 2017 the Vize T am is continously below -5 • C after 5 November. It is never above 0 • C. At the 115 Popova station T am very close to or above 0 • C during first two weeks of October. It is mostly below -5 • C after 7 November.
In the Barents Sea at the Kongsøya stations T am is longer periods below -5 • C only in 20-30 November, and 25 Dec onwards.
At the Varandey station T am is above 0 • C in 1-10 October, and continously above -5 • C up to 7 November. Afterwards there are periods up to 10 days with T am either above or below -5 • C. It seems that in the Kara Sea the time period after first week of November 2017 can be assumed to represent mostly cold winter conditions. In the Barents and Pechora Seas in November 120 and Dec there are also periods of cold winter conditions, but sometimes T am is between -5 and 0 • C, when there could have been moist snow effects on the radar signatures.

Data
In this section we present the data sets used in our study. CS-2's primary payload is the Synthetic Aperture Interferometric Radar Altimeter operating in the Ku-Band (13.6 GHz) which allows a much smaller sampling footprint (about 250 m in the satellite track direction) than earlier satellite radar altimeters.
Over sea ice, CS-2 echoes are assumed to scatter from the interface between the ice surface and the layer of overlying snow (Laxon et al. , 2013), thus enabling freeboard estimation, and further, sea ice thickness calculation by known snow depth and density. 130 The quantity altimeters measure is surface elevation, from which freeboard (F ) can be derived from. Assuming hydrostatic equilibrium and known sea ice density (ρ i ), snow density (ρ s ) and water density (ρ w ) the following equation can be formed: As can be seen from equation above, snow estimate has a large effect on retrieved SIT. In the case of radar altimeters, such as CS-2 used here, h s is also required for signal propagation speed correction in deriving F . A more detailed analysis on the 135 effect of snow can be found in (Kern et al., 2015). In this study, like in many other CS-2-based sea ice thickness products, climatological snow depth and density based on the Warren climatology (Warren et al. , 1999) are used.
In this study we have computed sea ice thickness with FMI implementation of the python sea-ice radar altimetry toolbox (PySiral, available in: https://github.com/shendric/pysiral). The primary input has been the European Space Agency's CS-2 Baseline-C Level-1B product. Our processing follows the algorithm of (

Reference Data
In our study we use the Russian AARI ice charts, ORAS5, PIOMAS, and TOPAZ4 ocean-sea ice reanalyses SIT data, and 150 MODIS daily SIT charts as reference data. They are used to evaluate the performance of our CS2-SAR SIT estimation method.
The TOPAZ4 SIT model data are also used to remap the CryoSat-2 training data set based on cumulative SIT distributions.

AARI ice charts
Arctic and Antarctic Research Institute (AARI) produces weekly ice charts for the Barents and Kara Seas (Bushuev et al. , 2007). The ice charts are based on the available visible, thermal infrared and SAR imagery and reports from coastal stations 155 and ships. The polygonalization of images and subsequent interpretation and mapping of ice conditions are carried out by skilled ice analysts. The main purpose of the weekly ice chart is to show the spatial distribution and characteristics of sea ice.
Many sea ice details within polygons are ignored. The AARI ice charts are in the digital SIGRID-3 vector file format (JCOMM , 2014b). We have not found any publication that discusses the uncertainty of the AARI weekly ice chart. The AARI ice charts are available at: http://wdc.aari.ru/datasets/d0004/.

160
The ice charts convey their information in codes that are explained in (JCOMM , 2014a reference. Looking at a satellite image, ice expert matches brightness of drifting ice outside the shore to brightness of land-fast ice with known thickness." Thus, we interpret the AARI ice thickness as the thickness of level ice in the area -analogous to the thermodynamic grown fast ice. For validation, polygons from AARI charts were gridded to 1 km by 1 km grid covering our study area. Then for each AARI ice chart polygon a mid-range SIC and mid-range SIT or upper or lower SIT values of the first, second, and third thickest ice 175 types were assigned. However, for FYI general type a SIT value of 95 cm is used (same as for thick FYI). An average SIT for each chart polygon was calculated as a sum of the SIT values for the three ice types weighted with their concentrations. Finally, the gridded charts for the Barents and Kara Seas were combined together, see an example in Fig. 8(c). This processing of the AARI charts has been used previously by . The total number of the AARI weekly charts used in the validation was 37 (time span from 4 Jan until 26 Dec, excluding the summer period) for both the Barents and Kara Seas. TOPAZ4 is a coupled ocean-sea ice model with advanced data assimilation system for the North Atlantic Ocean and Arctic (Sakov et al. , 2012). The data assimilation is based on the use of ensemble Kalman filter (Evensen , 1994 (Ricker et al. , 2017), based on fusion of Cryosat-2 and SMOS SIT is assimilated to the TOPAZ4 reanalysis SIT product.
The modelled sea-ice thickness provided by the TOPAZ4 CMEMS product is a diagnostic variable based on prognostically calculated sea-ice volume divided by sea-ice concentration for each model grid cell. Accordingly, it provides the mean sea-ice 190 thickness for the portion of grid cell that is covered by sea ice.

ORAS5
Because TOPAZ4 seems to underestimate the Arctic ice thickness ( the thinnest ice CS2 is able to reliably measure. As TOPAZ4, the ORAS5 product provides the diagnostic sea-ice thickness for each grid cell that has been calculated by dividing the modelled sea-ice volume by concentration.

PIOMAS
Additionally we utilized the PIOMAS model (Zhang and Rothrock , 2003) SIT data which have a coarser resolution than comparisons, but we utilise it in defining the mapping of CS2 thicknesses to model compliant thicknesses in Section 3.

MODIS Daily SIT Chart
Daily MODIS ice thickness (h iM ) charts for our study area have been processed for two winters, from November 2015 to April 2017. All the details on the daily chart processing can be found in . The daily charts are based on all available Aqua and Terra MODIS h iM swath charts. Charts for very cloudy days were manually excluded. The The MODIS based SIT charts have pixel size of 1 km, but the cloud mask has 10 km pixel size. The daily h iM chart shows daily median h iM (h d iM ) for pixels which had at least two h iM samples from the swath charts. The requirement for having at least two valid h iM samples decreases the errors due to the misdetected clouds in the swath charts.

S-1 preprocessing
We georectified and sampled the S-1 SAR data into a 100-m pixel size. Absolute calibration was applied to get the backscat- DN i is the original image pixel value at location i, n i is the provided noise data at location i,and A i the provided calibration 235 coefficient at location i.
The S1 data was preprocessed by applying a linear incidence angle correction  to the HH channel and a combined incidence angle and noise floor correction to the HV channel, for details, see (Karvonen , 2017). The HH channel incidence angle correction has been tuned for sea ice. This leads to reduced performance on open water because of the varying contribution of the water surface due to waves (Bragg scattering) to the backscattering measured by SAR. After 240 incidence angle and noise floor corrections the image data were geo-rectified into the polar stereographic projection specified in Section 2. After geo-rectification the imagery were down-sampled to 500 m resolution, and finally the daily mosaics were constructed by overlaying the newer images over the older ones such that at each mosaic grid cell (pixel) the newest SAR data prior to the mosaic time label, which was defined to be 12:00 UTC daily here, was available.
The mosaic was initialized only in the beginning of the mosaicking (in this case in the beginning January 2016). In practice 245 the data at a given grid cell location were never older than three days from the mosaic time label. Separate mosaics for HH and HV channels were constructed. A land mask based on the GSHHG coastline (Wessel and Smith , 1996) was applied to the mosaics to exclude land areas. The SAR mosaics of 21 February 2017 for HH and HV channels are shown in Fig. 4 and a detail of the channel images in Fig. 5.
SAR images were processed to 8 bpp images by scaling the σ 0 between 1 and 255 (0 representing background). The scaling 250 for the HH channel is such that -30 dB corresponds to the pixel value of one and 0 dB corresponds to the pixel value of 255. For HV channel the decibel values are -40 and 0, respectively. For segmentation the meanshift algorithm (Fukunaga and Hostetler , 1975) was first applied to locate the modes of the two-dimensional (HH and HV) SAR data. The meanshift algorithm has been empirically adjusted so that about 10-15 modes will be produced after convergence. The initial 10-15 categories based on MS were then used as a starting point for iterated conditional modes (ICM) segmentation (Besag , 1986). A block diagram of the algorithm is presented in Fig. 6. In the first phase the CS-2 SIT values are assigned to SAR segments of the segmented SAR mosaic. In the next phase the difference function is computed for all the segments without an assigned SIT value yet between all the segments with an already assigned SIT value within range. In the following phase a SIT is assigned 265 to the segments without a SIT value using a median of SIT values of M segments with the smallest difference function values.
As a last step, a remapping based on histograms of the estimated SIT and PIOMAS, ORAS5 and TOPAZ4 model reanalyses SITs for the training data are performed to get a bias-corrected SIT estimate for each segment.
The SAR features were computed within a round-shaped window with a radius R. In this study we have used R=5 pixels.
Median of the feature values within each segment were then assigned to the segments as segment-wise texture features. The 270 features used in this study are: • HH local autocorrelations (C A ) • HV local autocorrelations (C A ) • HH entropies (E) • HV entropies (E) • HH local variogram slopes (V 1 )

280
• HV local variogram slopes (V 1 ) • HH edge point densities (D E and D C ) • HV edge point densities (D E and D C ) • HH corner point densities (D E and D C ) • HV corner point densities (D E and D C )

285
• HH/HV channel cross-correlation (C c ) The feature computation windows were overlapping by half of the window size (i.e. window radius) in the both the coordinate directions.
Entropy E (Shannon , 1948) was computed as

290
where p k 's are the proportions of each gray tone k within each computation window. Auto-correlation C A (Box and Jenkins, 1976) was computed as where I(k, l) is the pixel value at image location (k,l). Mean over the horizontal, vertical and diagonal directions i.e. (k, l) = (0, 1), (k, l) = (1, 0), (k, l) = (1, 1) and (k, l) = (1, −1) was used to accomplish directional independence. The computation 295 window is here denoted by W. σ W and µ W are the mean and standard deviation within the window, respectively.
The cross-correlation C c (Knapp and Carter , 1976) between the SAR polarization channels (HH and HV,) here denoted by X (HH) and Y (HV), is

300
where k and l refer to the row and column coordinates of the window centre image pixel, respectively, σ x and µ x are the mean and standard deviation, respectively, of the window in X and σ y and µ y are the mean and standard deviation, respectively, of the window in Y . N p is the number of pixels within the window denoted by W . We also computed texture features based on local variograms. The variograms were locally estimated in a window with a radius of five pixels. Assuming a stationary and isotropic process, the variogram γ is (locally) dependent on the inter-distance, 305 here denoted by h, only (Cressie , 1993), and can be estimated as where z i and z j are the pixel values at locations i and j, whose distance is h. We have computed the length of approximately linear part of the variogram as a function of h, and the slope of a linear fit of this linear part of the variogram. These are referred here as features V ch 1 and V ch 2 , where ch (channel) is either HH or HV.

310
The coefficient of variation is computed as where σ and µ are the standard deviation and mean within the window w. Based on the segmentation result and the SAR texture features complemeted by σ 0 HH and σ 0 HV the segment-wise median values of each texture feature and channel-wise SAR backscatterting coefficients were calculated, resulting to a total of 15 320 SAR features (13 texture features and two backscattering coefficients) for each segment. As an example of the segmentation the segmentation result of the SAR mosaic of 21 February 2017 with σ 0 HH and σ 0 HV medians assigned to SAR segments is shown in Fig. 3 and a detail of these HH and HV channel figures in Fig. 4. We use the SIT estimation and reference SIT data of this day as an example in the following sections.
According to our analysis the SIT difference between two segments was mainly explained by L1 difference of a few features 325 (σ 0 HV , HV entropy and HH edge density) and the spatial distance between segment means, based on a least squares fit of the training data (Jan-Apr and Oct-Dec 2016), but other features also had minor effect. For this reason we have used all the abovementioned features in our difference function because the number of features did not produce any computational problems with respect to execution time or hardware resources on a common desktop personal computer.
The 2016 (January-April, October-December) CS-2 and SAR data were used for training and 2017 (January-April, October-330 December) data for testing the method. The segment difference function T was defined as linear combination of the SAR feature differences and temporal and spatial distance between as segment pair: ∆t is the (absolute) time difference in days, ∆D is the distance difference between the centres of the segments, ∆f k is the difference between the SAR feature f k (k = 1, ..., 15) of the two segments in comparison. 2016 CS-2 thickness and SAR 335 mosaics were used in defining the coefficients c t , c d , c k , using a non-negative least squares (LS) fit. Non-negative LS was used because all the absolute differences should have an increasing effect on T .
As training data we used the CS-2 SIT values assigned to the segments. Each 2016 training data SAR segment with an assigned SIT were included and the neighboring assigned SIT values and the corresponding segment pair-wise feature differences were used in the LS fit, i.e. each segment pair with an assigned CS-2 SIT were utilized.
where sup x is the supremum within the range of x (SIT in our case) and F t (x) and We have compared the SIT values produced by the algorithm to the SIT based on the AARI ice charts, ORAS5, PIOMAS and TOPAZ4 reanalysis SITs to evaluate the performance of the SIT estimation. In the comparisons we have used the following measures of difference between the SIT estimates and the reference SIT:

385
The correlations for the CS-2 SIT assigned to SAR segments was somewhat higher, around 0.35. The correlations w.r.t. the two reference data sets for each 2017 winter month and their averages and standard deviations are presented in Table I, the corresponding L1 differences in Table II, the biases in Table III and the RMSD's in Table IV. Also the differences and correlations between the three reference SIT data sets are shown in Table V. CS-2 in the tables refers to the CS-2 data mapped to the segments with more than N CS-2 measurements available. The biases (i.e. signed L1 differences) in Table III  It should be noted that also the reference SIT data sets differ from each other. These differences are shown in Table V Furthermore, the reader should note that in our evaluation the compared randomly sampled points were located in highly dissimilar regions, characteristic for the different data sets, with different sizes and shapes (such as ice model grid points, SAR segments, ice chart polygons).
In order to assess how far the SIT can be interpolated and still give usable estimates, we studied the effect of distance and 435 time on the SIT difference (or in other words, estimation error) using our training data set of 2016. We used the assigned CS-2 SIT values as reference and searched for sets of segments with either constant time difference or constant distance and defined the increase of ice thickness difference as a function of the difference in the other.
The average increase in estimation error for the 2016 training data as a function of time was 1.6 cm/day and of space was 1.3 cm/100 km. This indicates that the contribution of distance difference to the total difference is in maximum (corresponding 440 to the search range boundaries) around 13 cm for the 1000 pixel spatial search range and 16 cm for the 10 day temporal search range. An example of the SIT estimation of 21 February 2017 without and with the remapping can be seen in Fig. 8. For reference also the AARI SIT and the TOPAZ4 and ORAS5 reanalysis SIT have been included in Fig. 8. This figure represents a typical case with a general agreement of the thinner and thicker ice fields generally in agreement compared to the reference data but still indicating quite much difference due to the details of the local SIT distribution, given in higher level of detail in 445 the estimated CS2/SAR SIT.
As we are using a (segment) difference function, we compute its median value for each segment mapped. This value corresponds to the mapping and the larger it is the larger the uncertainty of the estimation can be considered. We have scaled the difference to the range 0-100 (based on the training data difference function values) and in Fig. 10 the scaled difference function for the estimation of 21 February 2017 is shown. The highest uncertainty is found at the ice edge and in the river 450 Ob delta area. The uncertainly could probably also be used for iteratively re-estimating the SIT of segments with the highest uncertainties.
We also performed some comparisons of the CS-2/S-1 SIT estimates against the daily MODIS SIT charts. The MODIS product only gives thin ice thicknesses and the thicker ice thicknesses are very roughly categorized, i.e. one category for SIT of 31-50 cm and another for SIT over 50 cm. There were MODIS SIT charts for the winter 2017 (January-April). The monthly 455 correlations between the thin MODIS SIT (1-30 cm) and the SIT estimated by the proposed method varied between 0.15 to 0.25. We also computed the averages of the estimated SIT for the three MODIS SIT categories, i.e. 1-30 cm, 31-50 cm and over 50 cm for each month and for the estimated SIT and estimated remapped SIT. These are presented in Table 5. It can be seen that the averages increase towards a thicker MODIS thickness category for both the estimates, but the increases are quite small and the averages are too high for both the methods, some less for the remapped estimate. It can also be seen that the averages 460 increase as a function of time, i.e. as ice gets thicker in the course of time. Based on this experiment we can conclude that the estimates, even with the proposed remapping, tends to overestimate specifically thin ice thickness. For visual evaluation and comparison to the SIT estimates of Fig. 9 we also show the MODIS SIT collage for February 21 2017 in Fig. 10. The color scale is the same as in Fig. 9, except that the thickest MODIS SIT category (SIT>50cm) is cyan in the figure. Because in a daily MODIS SIT chart SIT can typically be computed only for small areas due to cloud cover, the collage is composed of two In this study an algorithm for inter/extrapolating of the CS-2 SIT data over the S-1 daily SAR mosaic was developed and evaluated by comparisons to SIT derived from Russian Arctic AARI ice charts and two ice model reanalysis data sets. We wanted to demonstrate the potential of our method as well as to evaluate its performance. We found significant differences (both low correlation and high bias) between non-interpolated CS-2 SIT and the reference data. The match between CS2 and reference data improved with the introduction of segmentwise medians. However, there was still a significant difference, 480 especially a high positive bias due to different nature of CS-2 estimates and the reference data sets. To reduce this bias, we performed a mapping based on the PDF's (histograms) of the CS-2 SIT and PIOMAS/ORAS5 model reanalysis SIT of the training data set. This mapping reduced the positive bias. However, the remapped CS-2 data should not be understood as more correct, but as one more akin to model SIT.
An interesting result was that the inter/extrapolation does not decrease the correlation compared to (incomplete) SIT result-485 ing from CS-2 data assigned to SAR segments, or CS-2 data as individual measurements. In many cases correlations were even slightly higher for the results with the SAR data and inter/extrapolation included. Here the mapping was based on PI-OMAS/ORAS5 model reanalysis data of the training data of 2016. TOPAZ4 reanalysis still seems to underestimate the high SIT values, this is suggested e.g. by the comparison results between TOPAZ4 and ORAS5 SIT in Table V. In general, evaluation of ice thickness over Arctic is difficult because of lack of reliable reference data. Another problem is 490 the scale of the SIT measurements, estimates and models. SIT measurements, if available, are typically point measurements or measurements with a footprint of a few meters to a few tens of meters (e.g. SIT measurements by drilling or Electromagnetic induction based EM measurents), and the model grid cells, EO measurements or segment-wise SAR estimates represent averages over several sqare kilometers (km 2 ). Also ice charts give averages of typically quite large polygons of several km 2 .
We selected parameters, such as N, empirically. The sensitivity of SIT to parameters chosen was not studied thoroughly yet. assigned to them to get inter/extrapolated SIT for all the segments of a daily SAR mosaic. Currently the value of N is relatively small for a typical SAR segment. A higher value of N would be preferable. However, a higher value of N in turn would reduce the number of segments included in the estimation and also increase the possibility of too few assigned SIT values for the inter/extrapolation step. For this reason we used a moderate value of N.

500
The relatively larger differences w.r.t. the reference data sets can at least partly be explained by the different resolutions and level of detail of the SIT estimates. The CS-2/S-1 SIT has a resolution of 500 m (the SAR mosaic resolution) which is significantly higher than the resolution of the ORAS5, PIOMAS and TOPAZ4 reanalyses (over 10 km) and also the ice chart In the future we plan to study segment clustering and then assigning CS-2 SIT value median or mode of each segment 510 cluster to all segments of a segment cluster instead of single segment medians. An easy solution would be to merge small segments into their neighboring larger segments to remove all the segments smaller than a given area. This would give more CS-2 measurements within the segments but on the other hand also giving SIT estimates representing averages of larger areas.
Despite of this it should be noted that the segment boundaries, representing different ice fields, are still represented in the resolution of the segmentation. With these approaches we aim to overcome the problem of either having too few CS-2 SIT 515 values assigned to segments or too few segments with an assigned SIT. Also more detailed utilization of SIT distributions of segment clusters with a large enough number of SIT samples for forming statistically reliable SIT distributions should be studied. Including thin ice detection (thickness < 20 cm) from microwave radiometer using a method in (Makynen and Simila , 2019) could also improve the SIT estimation results. In some cases also large SAR segments could be split into smaller ones to get more detailed and more accurate local SIT estimates.

520
To match segments we compute the value of the segment-wise difference function T for the matched segments. The value of T could possibly be used to implement an iterative SIT assigning scheme. In this scheme only SIT values with T less than a selected threshold value would be assigned to segments and the process iterated using all the SIT values assigned during the previous iterations in assigning SIT values to the segments still without an assigned SIT at each iteration. The threshold to be applied could also be varying as a function of the iteration.

525
The evaluation of this study was performed for January-April and November-December. This was because both the CS-2 SIT and AARI ice chart ice stage of development were not available for wet surface conditions. For such conditions reliable estimation of SIT is difficult or even impossible with the current data and algorithms. This is due to the deviating behaviour of radar signal at the wet air/snow-ice boundary, compared to a frozen surface.
It may be possible to use the dependence of the biases on the phase of the winter (related to average SIT), such as early 530 freeze-up, freeze-up, melt-donw phase, to reduce the bias by applying a bias-reducing mapping dependent on the phase of the winter instead on a common mapping for the whole winter. The phase of the winter could roughly be estimated based on the estimated average ice thickness.
We also studied use of regularized linear regression (LASSO) (Tibshiran , 1996) for reducing the number of texture features needed in the SIT estimation. According to our first tests the performance was nearly similar with less texture parameters 535 involved, but on the other hand the need for computational resources with the amount of LS fit parameters used now were not significantly larger than with a reduced number of parameters, and thus we here used all the studied parameters in the LS fit in this study. In the future reduction of LS fit parameters could be studied more if necessary e.g. for more efficient computing to cover larger sea areas in a reasonable time. This may be required to meet conditions for operational production, e.g. for a pan-Arctic high-resolution SIT product.

540
Our algorithm can be easily adapted to any satellite altimeter SIT product. It is especially relevant for the future CRISTAL mission (Kern et al., 2020). After its expected launch in late 2027, CRISTAL shall provide a time critical SIT product, which can be merged with SAR data to interpolate SIT between the CRISTAL ground tracks. Thus our study is one of the necessary steps in introducing satellite altimeter data into the realm of operational ice charting. Before CRISTAL, our algorithm can be applied to ICESat-2 as well as to Sentinel-3 based SIT estimates.