3.1 CryoSat-2 data

CS2’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 300 m in the satellite along-track direction; Scagliola, 2013) than traditional pulse-limited altimeters. Over sea ice, CS2 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, denoted by F here, estimation, and further SIT calculation assuming known snow depth and density and sea ice density (e.g. climatologically derived).

The quantity altimeters measure is surface elevation, from which freeboard can be derived. Assuming hydrostatic equilibrium and known sea ice density (ρ_i), snow density (ρ_s) and thickness (h_s), and water density (ρ_w) the following equation can be formed for the SIT (h_i) estimation:

$$\begin{array}{}\text{(1)}& h_{\mathrm{i}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{i}}}}F-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{i}}}}{h}_{\mathrm{s}}.\end{array}$$

In dry snow conditions, the radar scattering and backscattering penetrate the snow cover, and the backscattering comes from the ice surface; i.e. the ice freeboard is measured by a radar altimeter. In dry snow conditions h_i can be estimated from the ice freeboard F_i:

$$\begin{array}{}\text{(2)}& h_{\mathrm{i}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{i}}}}{F}_{\mathrm{i}}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{i}}}}{h}_{\mathrm{s}}.\end{array}$$

As can be seen from equations above, the h_s estimate has a large effect on retrieved SIT. In the case of radar altimeters, such as CS2 used here, h_s is also required for wave propagation speed correction in deriving F. A more detailed analysis on the effect of snow can be found in Kern et al. (2015). In this study, like in many other CS2-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 SIT with FMI implementation of the Python sea ice radar altimetry toolbox (pysiral available in https://github.com/shendric/pysiral, last access: 10 May 2022). The primary input has been the European Space Agency's CS2 Baseline-D Level-1B product. The processing follows the algorithm of Ricker et al. (2014) and Hendricks et al. (2021a). For auxiliary data we have used the DTU15 mean sea surface height product, EUMETSAT OSISAF sea ice concentration and sea ice type, and Warren et al. (1999) snow depth and density data with the 50 % reduction over first-year ice first proposed by Kurtz and Farrell (2011).

3.2 Sentinel-1 SAR

All the available Copernicus S-1 C-band dual-polarized extra wide (EW) swath mode level 1 ground-range-detected medium-resolution (GRDM) SAR data with the HH–HV polarization channels over the BKS study area and period (January–April and October–December 2016 and January–April and October–December 2017) were used in this study. The S-1 SAR data are publicly available through the European Space Agency (ESA) Copernicus Science Hub (https://scihub.copernicus.eu/, last access: 10 May 2022). The S-1 EW SAR images cover a region of 410 km by 400 km in size, with a pixel size of 40 m and a spatial resolution of 95–91 m by 90 m (range by azimuth), and they have θ₀ variation from 19 to 47^∘.

3.3 Reference data

In our study we use the Russian Arctic–Antarctic Research Institute (AARI) ice charts, Ocean Reanalysis System 5 (ORAS5), Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS), and TOPAZ4 ice model ocean–sea ice reanalyses SIT data, CS2SMOS SIT data (Ricker et al., 2017), and MODIS daily SIT charts (Makynen and Karvonen, 2017a) 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.

3.4 S-1 preprocessing

We georectified and sampled the S-1 SAR data into a 100 m pixel size. The calibration to provide the backscattering coefficients (${\mathit{\sigma }}_{\mathrm{HH}}^{\mathrm{0}}$ and ${\mathit{\sigma }}_{\mathrm{HV}}^{\mathrm{0}}$ for HH and HV channels) at the image location (pixel i) was performed according to the following equation (Bourbigot et al., 2016):

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_i^{\mathrm{0}}={\displaystyle \frac{\left({\mathrm{DN}}_i^{\mathrm{2}}-n_i\right)}{A_i^{\mathrm{2}}}}.\end{array}$$

DN_i is the original image pixel value at location i, n_i is the provided noise data at location i, and A_i is the provided calibration coefficient at location i.

The S1 data were preprocessed by applying a linear incidence angle correction (Makynen and Karvonen, 2017b) 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 backscatter by ocean waves. After incidence angle and noise floor corrections the image data were georectified into the polar stereographic projection specified in Sect. 2. After georectification the images 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, were available.

Figure 4Sentinel-1 HH (a) and HV (b) mosaic σ⁰ segment median on 21 February 2017.

Figure 5Cropped area of the HH (a) and HV (b) segment median images of Fig. 4.

The mosaic was initialized only in the beginning of the mosaicking (in this case in the beginning January 2016). In practice the data at a given grid cell location were never older than 3 d from the mosaic time label. Separate mosaics for the HH and HV channels were constructed. Finally, a land mask based on the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) coastline (Wessel and Smith, 1996) was applied to the mosaics. The SAR mosaics of 21 February 2017 for the HH and HV channels are shown in Fig. 4 and small parts of the mosaics (more details visible) in Fig. 5.

SAR images were processed to 8 bits-per-pixel images by scaling the σ⁰ between 1 and 255 (0 representing background). The scaling for the HH channel is such that −30 dB corresponds to the pixel value of 1 and 0 dB corresponds to the pixel value of 255. For the HV channel the decibel values are −40 and 0, respectively. For segmentation the meanshift (MS) 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 mode (ICM) segmentation (Besag, 1986).

We also identify the low-SIC areas (SAR segments with SIC < 50 %) and exclude them from the SIT estimation procedure. The algorithm utilizing both SAR and AMSR2 microwave radiometer data of Karvonen (2017) is used to locate the low-SIC areas.

4.1 Interpolation of CS2 SIT using S-1 SAR

In this subsection we describe the mapping and interpolation of the CS2 SIT to SAR segments. SAR segments represent uniform areas in the SAR image, and each SAR segment is here assumed to have a uniform SIT value derived from the CS2 SIT values within the segment or interpolated from the CS2 SIT values of other SAR segments with CS2 SIT values in case there are not enough CS2 SIT values within a SAR segment.

As we want to provide SIT estimates for uniform areas of the SAR mosaics, i.e. SAR segments, the mosaics first need to segmented. After segmentation the CS2 SIT values are mapped to SAR segments where enough single CS2 SIT measurements are available. The SIT value assigned to these segments is the median value of the CS2 SIT values within the segment. Then the SIT is estimated (interpolated) for the rest of the segments by utilizing a difference function describing the similarity of the segments and pairwise distance of each pair of segments in time and space. The difference function includes difference of several SAR texture features in addition to the temporal and spatial pairwise distance of the segments. The data of January–April and October–December 2016 were used for training the algorithm and January–April and October–December 2017 data were used for evaluating the algorithm performance.

The SAR features were computed within a round-shaped window with a radius R. In this study we have used R=5 pixels. The median of the feature values within each segment was then assigned to the segments as segmentwise texture features. The features used in this study are as follows.

• HH backscattering coefficient (${\mathit{\sigma }}_{\mathrm{HH}}^{\mathrm{0}}$ in decibels)

• HV backscattering coefficient (${\mathit{\sigma }}_{\mathrm{HV}}^{\mathrm{0}}$ in decibels)

• HH entropies (E_HH)

• HV entropies (E_HV)

• HH local autocorrelations ($C_{\mathrm{A}}^{\mathrm{HH}}$)

• HV local autocorrelations ($C_{\mathrm{A}}^{\mathrm{HV}}$)

• HH/HV channel cross-correlation (C_c)

• HH local variogram slopes ($V_{\mathrm{1}}^{\mathrm{HH}}$)

• HV local variogram slopes ($V_{\mathrm{1}}^{\mathrm{HV}}$)

• HH coefficient of variations ($C_{\mathrm{v}}^{\mathrm{HH}}$)

• HV coefficient of variations ($C_{\mathrm{v}}^{\mathrm{HV}}$)

• HH edge point densities ($D_{\mathrm{E}}^{\mathrm{HH}}$)

• HV edge point densities ($D_{\mathrm{E}}^{\mathrm{HV}}$)

• HH corner point densities ($D_{\mathrm{C}}^{\mathrm{HH}}$)

• HV corner point densities ($D_{\mathrm{C}}^{\mathrm{HV}}$)

The feature computation windows were overlapping by half of the window size (i.e. window radius) in both coordinate directions.

Entropy E (Shannon, 1948) was computed as

$$\begin{array}{}\text{(4)}& E=-\sum _{k=\mathrm{0}}^{\mathrm{255}}{p}_k{\log }_{\mathrm{2}}{p}_k,\end{array}$$

where p_k values are the proportions of each grey tone k within each computation window. Auto-correlation C_A (Box and Jenkins, 1976) was computed as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& C_{\mathrm{A}}(k,l)=\\ & {\displaystyle \frac{{\sum }_{ij\in W}\left(I(i-k,j-l)-{\mathit{\mu }}_W\right)\left(I(i,j)-{\mathit{\mu }}_W\right)}{\left|B\right|{\mathit{\sigma }}_W^{\mathrm{2}}}},\end{array}\end{array}$$

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)=(\mathrm{0},\mathrm{1})$, $(k,l)=(\mathrm{1},\mathrm{0})$, $(k,l)=(\mathrm{1},\mathrm{1})$ and $(k,l)=(\mathrm{1},-\mathrm{1})$, was used to accomplish directional independence. The computation 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

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& C_{\mathrm{c}}(k,l)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}}{\mathit{\sigma }}_x{\mathit{\sigma }}_y}}\\ & \times \sum _{i,j\in W}\left(X(k+i,l+j)-{\mathit{\mu }}_y\right)\left(Y(k+i,l+j)-{\mathit{\mu }}_y\right),\end{array}\end{array}$$

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 5 pixels. Assuming a stationary and isotropic process, the variogram γ is (locally) dependent on the inter-distance, here denoted by d, only (Cressie, 1993) and can be estimated as

$$\begin{array}{}\text{(7)}& \mathit{\gamma }\left(h\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\left|N_d\right|}}\sum _{i,j\in N_d}|z_i-z_j{|}^{\mathrm{2}},\end{array}$$

where z_i and z_j are the pixel values at locations i and j, whose distance is d; N_d is the set of pixels with the mutual distance d; and $\left|N_d\right|$ is the number of these pixels within the window. We have computed the length of an approximately linear part of the variogram as a function of h and the slope of a linear fit of this linear part. These are referred to here as features $V_{\mathrm{1}}^{\mathrm{ch}}$ and $V_{\mathrm{2}}^{\mathrm{ch}}$, where ch (channel) is either HH or HV.

The coefficient of variation is computed as

$$\begin{array}{}\text{(8)}& C_{\mathrm{v}}={\mathit{\sigma }}_w/{\mathit{\mu }}_w,\end{array}$$

where σ and μ are the standard deviation and mean within the window w. C_v is computed separately for the HH and HV channels, respectively; i.e. we have $C_{\mathrm{v}}^{\mathrm{HH}}$ and $C_{\mathrm{v}}^{\mathrm{HV}}$.

Edge and corner points represent the locations with large sudden change in σ⁰, i.e. where high local gradients appear. Corner points are the points where the edge direction abruptly changes. The edge and corner point counts (N_e and N_c) were extracted for each segment using local binary patterns (Ojala et al., 1996) in a similar way as presented in Karvonen (2016). The proportion of the number of edge and corner points with respect to each segment area (in pixels) was computed for both polarization channels, and they were used as texture features. These features are denoted here by $N_{\mathrm{e}}^{\mathrm{HH}}$, $N_{\mathrm{e}}^{\mathrm{HV}}$, $N_{\mathrm{c}}^{\mathrm{HH}}$ and $N_{\mathrm{c}}^{\mathrm{HV}}$.

Based on the segmentation result and the SAR texture features complemented by ${\mathit{\sigma }}_{\mathrm{HH}}^{\mathrm{0}}$ and ${\mathit{\sigma }}_{\mathrm{HV}}^{\mathrm{0}}$, the segmentwise median values of each texture feature and channelwise SAR backscattering coefficients were calculated, resulting in a total of 15 SAR features (13 texture features and two backscattering coefficients) for each segment. As an example, segmented SAR mosaics with ${\mathit{\sigma }}_{\mathrm{HH}}^{\mathrm{0}}$ and ${\mathit{\sigma }}_{\mathrm{HV}}^{\mathrm{0}}$ medians assigned to SAR segments on 21 February 2017 are shown in Fig. 4, and details of these HH and HV mosaics are shown in Fig. 5. We use the SIT estimation and reference SIT data of this day as an example in the following sections.

The segment difference function T was defined as the linear combination of the SAR feature differences and temporal and spatial distance between a segment pair:

$$\begin{array}{}\text{(9)}& T=c_t\mathrm{\Delta }t+c_d\mathrm{\Delta }D+\sum _{k=\mathrm{1}}^{N_f}{c}_k\left|\mathrm{\Delta }{f}_k\right|.\end{array}$$

Δt is the (absolute) time difference in days, ΔD is the distance difference between the centres of the segments and Δf_k is the difference between a SAR feature f_k (k=1, …, 15) in the two segments. The 2016 CS2 thickness and SAR mosaics (training data) were used in defining the coefficients c_t, c_d and 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.

According to our analysis the SIT difference between two segments was mainly explained by the absolute difference, i.e. L1 difference, of a few features. ${\mathit{\sigma }}_{\mathrm{HV}}^{\mathrm{0}}$, HV entropy and HH edge density and the spatial distance between segment means were the most significant features based on a least-squares fit of the training data. However, other features also had a minor effect (small but non-zero coefficients based on the LS). For this reason we have used all the above-mentioned 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. Because the training data set consisted of data of a whole year with CS2 SIT data (January–April, October–December), we consider this to be a representative training data set and assume there is no significant overfitting.

As training data we used the CS2 SIT values assigned to the segments. Each 2016 training data SAR segment with an assigned SIT was included, and the neighbouring assigned SIT values and the corresponding segment pairwise feature differences were used in the LS fit; i.e. each segment pair with an assigned CS2 SIT was utilized.

SIT for a segment (S) without a SIT value was obtained as follows: T between S and all the other segments within a predefined time and distance range was computed, and the segments were ordered by ascending T. Then the CS2 SIT values of the other segments, in addition to possible CS2 SIT values with segment S (less than N values), were included in a set of SIT values, until the number of the values in the list exceeds the threshold value N (strictly equal to N or more than N). The SIT assigned to S was the weighted average of the SIT values of the collected SIT value set. The weights used in the averaging are relative to the inverse values of T corresponding to the segment from which the SIT values come. The possible (less than N) CS2 SIT values of the segment itself are included with the same weight as the CS2 SIT values within the segment with the smallest non-zero value of T. We have used parameter values of N=7 in this study.

Figure 7(a) Dependence of the SIT difference on the pairwise SAR segment texture difference, (b) dependence of the SIT difference on the pairwise SAR segment distance and (c) dependence of the SIT difference on the pairwise SAR segment time difference. The cyan lines in the figures represent the linear least-squares regression.

Download

In Fig. 7 the dependence of the SIT difference on the SAR feature linear combination, pairwise segment distance and pairwise segment time difference is shown for the 2016 data. Here the time difference is restricted to be less than 1 d (i.e. using only the same day data) and distance less than 500 pixels (250 km) for Fig. 7a, describing the dependence on the pairwise segment texture difference. For Fig. 7b, describing the dependence of the SIT difference on the pairwise segment distance, and Fig. 7c, describing the dependence of the SIT difference on the segment time difference, the normalized pairwise texture difference is restricted to be under 0.1. For Fig. 7b the time difference less than 1 d, and for Fig. 7c the distance difference is restricted to be less than 250 km. It can be seen that there is a clear correlation between the linear combination of the texture feature differences and the SIT difference. There is quite a large deviation; however, the correlation is around 0.5. For the pairwise segment distance and the pairwise segment time difference, correlations are smaller (0.28 and 0.22, respectively), but there is still a visible trend. Based on the slopes of the linear fits, the average increase in SIT difference is 2.9 cm (100 km)⁻¹ and 2.1 cm d⁻¹.

4.2 Ice chart and model-compliant CS2 product

After defining the coefficients c_t, c_d and c_k in Eq. (9) using the training data, we tested the SIT estimation for the 2016 training data. It was observed that SIT is significantly overestimated in the 2016 training data when compared to the modelled and AARI SIT. Therefore, we introduce a mapping based on the PIOMAS model reanalysis SIT data to reduce this SIT overestimation. The CS2 SIT measurements are now further mapped based on the training data and the corresponding PIOMAS reanalysis SIT. The coarser-resolution PIOMAS data were selected for this purpose because ORAS5 does not give reasonable SIT values for SIT less than 50 cm (Tietsche et al., 2017). The adapted SIT histogram mapping is based on the normalized histograms and minimizing the Kolmogorov–Smirnov distance (i.e. mapping the cumulative histograms) of the CS2/S-1 SIT and the TOPAZ4 reanalysis SIT. The Kolmogorov–Smirnov distance, D_m, for a variable x is defined as

$$\begin{array}{}\text{(10)}& D_m={\sup }_x\left|F_{\mathrm{t}}\right(x)-F_{\mathrm{c}}(x\left)\right|,\end{array}$$

where sup_x is the supremum within the range of x (SIT in our case), and F_t(x) and F_c(x) are the cumulative probability density functions (CDFs) of the PIOMAS SIT and CS2 SIT computed for the 2016 training data set. In our case we have quantized the SIT into integer centimetres with 16 bits, and the SIT histograms are used as approximations of probability density functions (PDFs).

This mapping based on the 2016 training data CS2 SIT and PIOMAS model SIT is shown in Fig. 8 as a green curve, For reference the AARI ice chart SIT, ORAS5 model SIT and TOPAZ4 model SIT (blue curve) matching results are also shown in Fig. 8. It can be seen that TOPAZ4 reanalysis SIT mapping is at a significantly lower level than the mapping for the other two models.

Figure 8The correspondence of the CryoSat-2 SIT to AARI ice chart SIT and modelled SIT from the three models studied based on histogram matching.

Download

5.1 Measures used in evaluation of the results

We have compared the SIT values produced by the CS2/S-1 algorithm to the SIT based on the AARI ice charts, ORAS5, PIOMAS and TOPAZ4 reanalysis to evaluate the algorithm performance. In the comparisons we have used the following measures of difference between the CS2/S-1 SIT estimates and the reference SIT:

$$\begin{array}{}\text{(11)}& {\displaystyle }C={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{s}}\mathit{\sigma }{\mathit{\sigma }}_{\mathrm{ref}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{s}}}\left(\left(X_i-\mathit{\mu }\right)\left(X_i^{\mathrm{ref}}-{\mathit{\mu }}_{\mathrm{ref}}\right)\right),\text{(12)}& {\displaystyle }{D}_{\mathrm{L}\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{s}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{s}}}\left|X_i-X_i^{\mathrm{ref}}\right|,\text{(13)}& {\displaystyle }{D}_{\mathrm{sgn}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{s}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{s}}}\left(X_i-X_i^{\mathrm{ref}}\right),\text{(14)}& {\displaystyle }{D}_{\mathrm{RMS}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{s}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{s}}}{\left(X_i^{\mathrm{est}}-X_i^{\mathrm{ref}}\right)}^{\mathrm{2}}}.\end{array}$$

N_s refers to the number of samples (number of grid points involved) used in the comparison, and X_i (i=1, …, N_s) represents the estimated values of SIT and $X_i^{\mathrm{ref}}$ represents the values of the reference SIT data at the same location as X_i. C is the correlation, D_L1 is the L1 difference and D_sgn is the signed L1 difference giving the estimation bias, with positive bias indicating overestimation and negative bias indicating underestimation. D_RMS is the root-mean-square difference. μ and μ_ref are the means of the estimated SIT and reference SIT, and σ and σ_ref are the standard deviations of the estimated SIT and reference SIT, respectively. All the difference measures were computed in the resolution of each reference data set; i.e. the CS2 SIT assigned to SAR segments and interpolated, given in a high-resolution grid, were down-sampled to each reference data set resolution. This approach was not applicable to the single CS2 SIT measurements, and for them the difference measures were computed directly, using the nearest reference data set grid point SIT for each measurement as the reference SIT.

As the CS2 SIT data cover only the time period from October to April, the weather conditions were mostly with dry snow conditions; i.e. T_a was mainly below 0^∘ for the time range studied. Only in the early winter (October–November) there were some periods with T_a above zero and the potential to have wet snow on sea ice. In this study we do evaluate the effect of wet snow on the SIT estimation. For reliable estimation of the effect of wet snow conditions on estimation, more data acquired during the melting period would be required.

5.2 Difference statistics between the estimates and reference data

The statistics from the comparisons between the CS2/S-1 SIT, without and with the proposed histogram mapping, and the different reference SIT data sets are shown in Tables 1–4. The correlations between different SIT data sets for each 2017 winter month and their averages and standard deviations are shown in Table 1, the corresponding biases in Table 2, the L1 differences in Table 3 and the RMSDs in Table 4. For reference, the average difference measures and correlations between the different reference SIT data sets are shown in Table 5.

Table 1Monthly average cross-correlations between the SIT estimates studied with respect to the reference data set SIT. SIT refers to the proposed algorithm SIT, RSIT refers to the proposed algorithm SIT with remapping based on the histogram mapping with model data, CS2 refers to the CS2 values assigned to the SAR segments (without inter- or extrapolation), A refers to the SIT derived from the AARI ice charts, T refers to the TOPAZ4 model reanalysis SIT, and O refers to the ORAS5 model SIT.

Download Print Version | Download XLSX

Table 2Monthly average bias in centimetres between the different SIT estimates with respect to the SIT reference data sets. The symbols are the same as in Table 1.

Download Print Version | Download XLSX

Table 3Monthly average L1 difference in centimetres between the different SIT estimates with respect to the reference SIT data sets. The symbols are the same as in Table 1.

Download Print Version | Download XLSX

Table 4Monthly average RMSD in centimetres between the different SIT estimates with respect to the reference SIT data sets. The symbols are the same as in Table 3.

Download Print Version | Download XLSX

Table 5Average difference measures between the AARI ice chart SIT (A), CS2SMOS SIT (CS2SM) and SIT of the three ice models. O refers to ORAS5, T refers to TOPAZ4 and PM refers to PIOMAS.

Download Print Version | Download XLSX

On average the correlation between the non-interpolated CS2 SIT estimates and the AARI ice chart SIT was very small, only 0.15. There was a significant bias, CS2 SIT being on average 83 cm larger than the AARI chart SIT. This also reflects the high L1 differences and RMSDs between the CS2 and AARI SITs. The statistics of the CS2 SIT with respect to the AARI SIT are given in the first columns of Tables 1–4. The numbers were similar with respect to the three ice models, showing small correlations and large positive bias (overestimation by CS2). Because of their similarity these values were not included to the tables.

The correlations for the interpolated CS2 SIT assigned to the SAR segments was significantly higher, around 0.64. The monthly correlations between the interpolated CS2/S-1 SIT and the remapped interpolated CS2/S-1 SIT with respect to the AARI SIT were similar and in the range 0.52–0.76. The highest values were reached in April and October. There CS2/S-1 SIT correlations were quite similar with respect to the model data for all the studied models.

The monthly signed L1 difference in Table 2 also increases for the estimated SIT with respect to both the AARI SIT and model reanalysis SIT as the ice gets thicker, indicating significant overestimation. For the remapped interpolated SIT, bias remains at a lower level for all the studied months. With respect to ORAS5 there was even small underestimation on average, but with respect to TOPAZ4 reanalysis SIT there was some overestimation. The averages and maximum biases for the interpolated SIT with respect to the AARI SIT were 27 and 46 cm (maximum reached in March). The corresponding values with respect to the AARI SIT for the remapped interpolated SIT were 6 and 14 cm (maximum underestimation in April). The corresponding values with respect to the TOPAZ4 reanalysis SIT were 29 and 49 cm and 8 and 16 cm, respectively. The corresponding values with respect to the ORAS5 SIT were 6 and 20 cm and −15 and −20 cm, i.e. slight underestimation. For the PIOMAS model data the corresponding figures were 14 and 31 cm and −7 and −13 cm.

The monthly L1 differences in Table 3 were smaller for the freeze-up period and higher for the winter months, starting to decrease in April. The L1 differences with respect to the AARI SIT and TOPAZ4, ORAS5 and PIOMAS reanalysis SIT were smaller for the remapped estimated SIT, which was expected as the target of the remapping was to reduce the relatively large positive bias. The average and maximum L1 differences with respect to AARI SIT were 32 and 51 cm for the interpolated SIT and 14 and 21 cm for the remapped interpolated SIT. The maximum values were reached in January and March. The corresponding values with respect to TOPAZ reanalysis were 32 and 52 cm and 17 and 26 cm, with respect to ORAS5 SIT 29 and 41 cm and 24 and 31 cm, and with respect to PIOMAS SIT 29 and 42 cm and 21 and 30 cm. The difference maxima were reached in February–April for all the cases.

The RMSDs in Table 4 have a similar behaviour as the L1 differences. RMSD also increases as the average ice thickness increases in the course of the winter, and the smallest RMSD values were reached during the freeze-up months (October–December). The average and monthly maximum RMSD values for the estimated SIT with respect to the AARI SIT were 56 and 81 cm, and the corresponding values for the remapped estimated SIT were 25 and 34 cm. The corresponding values with respect to the TOPAZ4 reanalysis SIT were 58 and 84 cm for the interpolated SIT and 29 and 41 cm for the remapped interpolated SIT. The corresponding values with respect to the ORAS5 SIT were 49 and 65 cm for the estimated SIT and 38 and 50 cm for the remapped interpolated SIT. The corresponding values with respect to PIOMAS SIT were the 52 and 73 cm for the estimated SIT and 35 and 46 cm for the remapped interpolated SIT. The average maxima were again reached during the winter months (January–April).

For the CS2 SIT the difference measures shown in Tables 2–4 indicate large positive biases (overestimation) and large L1D and RMSD with respect to the AARI SIT, TOPAZ4, ORAS5 and PIOMAS reanalysis SIT. These values are highest for the winter months with the thickest ice and smaller for the freeze-up/early winter (October–December). The cross-correlations with respect to all reference data sets were low, significantly lower than for the interpolated CS2/S-1 SIT.

It should be noted that the reference SIT data sets also differ from each other; see Table 5. The average cross-correlations between the reference data sets were in the range 0.54–0.60. The TOPAZ4 and AARI SIT values were on average quite close to each other, having only low bias, but the ORAS5 SIT was then on average 35–40 cm above them. Furthermore, it is noted 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 from the assigned CS2 measurement the SIT can be interpolated and still give usable estimates, we studied the effect of distance and time on the SIT difference (or in other words, estimation error) using our training data set. We used the assigned CS2 SIT values as reference and searched for sets of segments with either constant time difference or constant distance and defined the increase in SIT difference as a function of the difference in the other. The average increase in estimation error for the training data set as a function of time was 2.1 cm d⁻¹ and of distance was 2.9 cm (100 km)⁻¹. This indicates that the contribution of distance difference to the total difference is at a maximum (corresponding to the search range boundaries) around 14 cm for the 1000 pixel (500 km) spatial search range and 21 cm for the 10 d temporal search range.

Figure 9(a) CS2/S-1 SIT, (b) remapped CS2/S-1 SIT, (c) AARI ice chart SIT, (d) TOPAZ4 ice model reanalysis SIT, (e) ORAS5 model ice model SIT and (f) CS2/SMOS SIT for 21 February 2017.

An example of the SIT estimation on 21 February 2017 without and with the histogram remapping can be seen in Fig. 9. For reference the AARI SIT and TOPAZ4, ORAS5 reanalysis SIT, and CS2SMOS SIT have also been included in Fig. 9. This figure represents a typical case of the thinner and thicker ice fields generally in agreement with the reference data but still indicating differences due to the details of the local SIT distribution, given at a higher level of detail in the estimated CS2/S-1 SIT.

Table 6Comparison between the CS2 SIT assigned to SAR mosaic segments and CS2SMOS SIT.

Download Print Version | Download XLSX

Table 7Comparison between the CS2 assigned to SAR mosaic segments with the mapping and CS2SMOS SIT.

Download Print Version | Download XLSX

We also included a comparison of the proposed SIT estimates (interpolated and remapped versions) to the CS2SMOS SIT data (Ricker et al., 2017) for our 2017 test data set. These comparison results are shown in Tables 6 and 7. The results are discussed more in the following subsection.

5.3 Difference maps and estimation uncertainty

We also provide difference maps of the unmapped and remapped SIT with respect to the AARI ice chart SIT in Fig. 10 for the 21 February 2021 case and for the 2017 test data set on average. It can be seen that the unmapped interpolated SIT overestimates SIT with respect to the AARI ice chart SIT, for both the 21 February 2017 case and the average in many areas. For the remapped interpolated SIT estimates the overestimation is significantly smaller with respect to the AARI ice chart SIT for both the case of 21 February 2017 and on average.

Figure 10(a) Difference of the proposed (unmapped) SIT for 21 February 2017 and AARI ice chart SIT, (b) difference of the proposed remapped SIT and AARI ice chart SIT for 21 February 2017, (c) average difference of the proposed (unmapped) SIT and the AARI ice chart SIT for the 2017 data, and (d) average difference of the proposed remapped SIT and the AARI ice chart SIT for the 2017 data.

Figure 11(a) Difference of the proposed (unmapped) SIT and CS2SMOS SIT for the 21 February 2017 case, (b) difference of the proposed remapped SIT and CS2SMOS SIT for the 21 February 2017 case, (c) average difference of the proposed (unmapped) SIT and the CS2SMOS SIT for the 2017 data, and (d) average difference of the proposed remapped SIT and the CS2SMOS SIT for the 2017 data.

SIT difference charts with respect to the CS2SMOS SIT on 21 February 2017 and the average difference charts for 2017 are shown in Fig. 11. It can be seen that the SIT differences are quite similar to the differences with respect to the AARI ice charts. There is some more underestimation in the northeastern part of the study area compared to the difference maps with the AARI SIT. This comparison also indicates that the interpolated SIT estimates deem to overestimate SIT, and the remapped interpolated SIT values are less biased.

As we are using a (segment) difference function for each segment, this value corresponds to the mapping and the larger it is the larger the uncertainty of the SIT estimation can be considered. We have scaled the difference T values to the range 0–100 (based on the training data difference function values), and in Fig. 12 the scaled difference function for 21 February 2017 is shown. The highest uncertainty is found at the ice edge and in the Ob River delta area. The uncertainly could probably also be used for iteratively re-estimating the SIT of segments with the highest uncertainties.

Figure 12Segment difference (T) as a measure of uncertainty of the SIT estimates. This is the segment difference for the 21 February 2017 case. The values of T have been scaled to the range 0–100, based on the variation in T in the training data. The black area is open water.

5.4 Comparison to MODIS SIT

We also performed comparisons of the interpolated CS2/S-1 SIT estimates against the daily MODIS SIT charts. The MODIS product only gives thin-ice thicknesses, and the thicker ice thicknesses are only roughly categorized: SIT in the 1–30 cm range and then one category for SIT of 31–50 cm and another for SIT over 50 cm. The monthly correlations between the thin MODIS SIT (1–30 cm) and the SIT estimated by the proposed method varied between 0.15 and 0.25. We also computed the averages of the estimated SIT for the three MODIS SIT categories for each month and for the estimated SIT and estimated remapped SIT. These are presented in Table 8. 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 high for both the methods, some less for the remapped estimate. It can also be seen that the averages 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 SIT estimates, even with the proposed remapping, tend to overestimate SIT, specifically for thin ice.

Table 8Monthly averages of estimated SIT (SIT) and remapped SIT (Rem. SIT) for the MODIS daily SIT chart categories. All values are in centimetres. The MODIS SIT categories are 1–30 cm (Cat1), 31–50 cm (Cat2) and over 50 cm (Cat3).

Download Print Version | Download XLSX

For visual evaluation and comparison to the SIT estimates of Fig. 9, we also show the MODIS SIT collage for 21 February 2017 in Fig. 13. For better visual appearance the colour scale is different from that of Fig. 9. 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 2 weeks of MODIS SIT charts with the most recent available MODIS SIT value at each grid point. The areas with no data during the 2-week period are indicated by the black colour in the figure.When comparing the MODIS SIT collage of Fig. 10 with the daily CS2/S-1 SIT estimate in Fig. 8b, we notice that most of the ice-covered regions of the study area belong to the thickest MODIS SIT category (SIT > 50 cm). The thin-ice areas according to the MODIS SIT are located near the Kara Gate (southern Kara Sea) and near the ice edge in the northwestern part of the study area. These are in general in agreement with the CS2/S-1 SIT. Also, some thicker ice patches within or in the vicinity of the MODIS SIT thin-ice areas can be identified in the CS2/S-1 SIT. The general pattern detecting thinner ice in the south and near the ice edge is present in both CS2/S-1 and MODIS SIT. This general large-scale pattern can be seen throughout the winter. However, the SIT estimates include many anomalous local details due to the dynamic nature of the ice field and the fact that the data are not exactly simultaneous.

Figure 13MODIS SIT collage based on 2 weeks of MODIS daily SIT charts before 21 February 2017. The areas without data are indicated by the grey colour, and all the ice thicker than 50 cm is indicated by the green colour, corresponding to 60 cm in the colour map.