1. Introduction
Sea ice thickness (SIT) is an essential sea ice parameter useful for multiple purposes. Remote sensing of sea ice thickness involves measurement and estimation of SIT by various satellite and airborne sensors. Accurate and timely information on sea ice thickness is crucial e.g. for climate research, weather forecasting, navigation, and resource exploration in the polar regions (Ronkainen and others, Reference Ronkainen2018; Cheng and others, Reference Cheng2023). From sea ice thickness and sea ice concentration (SIC) also sea ice volume (SIV) can be derived in a straightforward manner. During the recent years, remote sensing techniques have played a significant role in enhancing understanding of sea ice properties and their spatial and temporal variations. In this introduction we first provide an overview of SIT estimation based or EO data to enlighten the context of this study. A comprehensive literature review of SIT estimation based on EO until 2017 has already been made for the Norwegian Defense Research Establishment in Hannevik (Reference Hannevik2017).
Satellite altimetry (SA) is one of the primary remote sensing techniques used to measure SIT. Altimetry satellites, such as CryoSat-2 radar altimeter (Laxon and others, Reference Laxon2013) and ICESat-2 laser altimeter (Kwok and others, Reference Kwok2019), utilize a radar altimeter (RA) or a laser altimeter (LA) to measure the ice surface height (RA) (Giles and others, Reference Giles, Laxon and Worby2008; Tilling and others, Reference Tilling, Ridout and Shepherd2018) or the joint height of the sea ice and the snow cover on sea ice (LA) (Kern and others, Reference Kern, Ozsoy-Çiçek and Worby2016; Kurtz and Markus, Reference Kurtz and Markus2012). Based on this information, and making assumptions of the snow cover thickness and the snow and ice densities, SIT can be estimated. SA instruments provide valuable data for mapping SIT over large areas. However, SA only provides measurements along narrow measurement lines. defined by the satellite track, and thus their spatial and temporal coverage is not especially suitable for areal high-resolution Near-Real-Time (NRT) SIT monitoring. In Kacimi and Kwok (Reference Kacimi and Kwok2022) ICESat-2 and Cryosat-2 data have been combined to get SIT measurements with a better spatial and temporal coverage compared to using data from an individual altimeter instrument.
Electromagnetic induction (EM) sensors are also used to measure sea ice thickness (Haas and others, Reference Haas, Lobach, Hendricks, Rabenstein and Pfaffling2009) of a measurement line along the instrument flight path. EM sensors for SIT measurements are typically mounted on an aircraft or a helicopter or outside of a ship sailing in ice. An EM sensor has two coils, a transmitter and a receiver coil, at each end. The primary EM field emitted by the transmitter coil mostly penetrates snow and ice and induces an eddy current near the sea-water surface under the ice. This is due to the high contrast in conductivity between resistive sea ice and conductive sea water. The secondary EM field caused by the eddy current is detected by the receiver coil. The ratio of intensities between the primary and secondary EM field is approximately related to the apparent conductivity. The apparent conductivity can then be transformed to the distance to the water surface. The distance from the instrument to the sea ice surface, or actually to the snow on the sea ice, can then be measured by a separate laser instrument, and the combined snow and sea ice layer thickness can then be derived from the EM distance and the distance to the snow/ice surface measured by the laser by subtraction. EM sensors provide high-resolution measurements but are limited to relatively small survey areas. EM data are typically used as reference data for EO observations or in limited-area scientific experiments.
Passive microwave radiometers (MWR's) can also be used to estimate SIT, reported e.g. in Martin and others (Reference Martin, Drucker, Kwok and Holt2004) and Tamura and Ohshima (Reference Tamura and Ohshima2011). Microwave radiometers measure the microwave radiation emission from the Earth's surface at specific microwave frequencies. Microwave frequencies (300 MHz–300 GHz) can penetrate sea ice to varying depths, allowing to estimate SIT based on the emitted microwave signal measured by the MWR. Microwave radiometers are particularly sensitive to the presence of liquid water within the ice. The brightness temperature measured by the radiometer is influenced by the emissivity of the ice, which is related to both the presence of liquid water and the thickness of the ice. Thicker ice generally has a higher brightness temperature compared to thinner ice. Depending on the frequency, the brightness temperature saturates at different ice thickness and estimation of thicker ice with thickness over the saturation thickness is not possible. The saturation ice thickness is larger for longer wavelengths. Satellite-based microwave radiometers, such as those onboard the Special Sensor Microwave Imager/Sounder (SSMIS) and the Advanced Microwave Scanning Radiometer for EOS (AMSR-E and AMSR2), have been used to estimate sea ice thickness. By combining MWR brightness temperature measurements at multiple MW frequencies MWR is capable of estimating SIT but only for relatively thin ice (Tamura and others, Reference Tamura2007; Nihashi and others, Reference Nihashi, Ohshima and Tamura2017; Demir and others, Reference Demir2022).
A passive microwave instrument proven to be useful for ice thickness estimation is the Soil Moisture and Ocean Salinity (SMOS) mission, developed by the European Space Agency (ESA). SMOS. is primarily focused on measuring soil moisture and ocean salinity. However, SMOS data can also be utilized to estimate SIT (Tian-Kunze and others, Reference Tian-Kunze2014; Gupta and others, Reference Gupta, Gabarro, Turiel, Portabella and Martinez2019). SMOS employs a passive radiometer that operates at L-band microwave frequencies (more accurately 1400–1427 MHz), which allows it to penetrate the ice surface layer and being able to measure the radiation emitted from the ice volume. By analyzing the brightness temperature data received from SMOS, information about the sea ice thickness can be retrieved. The estimation of sea ice thickness using SMOS data is based on the concept that thicker ice has a higher emissivity and therefore emits radiation differently compared to thinner ice. This difference in emitted radiation can be detected by SMOS and used to infer the thickness of the sea ice. Because of the limited penetration depth SMOS can only measure thickness of rather thin ice, typically less than 50 cm.
Combining data from both SMOS and CryoSat-2 radar altimeter can provide valuable synergistic information for estimating sea ice thickness. SMOS is particularly sensitive to thin ice and can provide information on areas with low sea ice thickness. SMOS data can also be used to identify regions where the ice cover is thin enough for the L-band microwave signals to penetrate through it and be detected. CryoSat-2, on the other hand, is designed to measure sea ice freeboard, i.e. the height of ice above the waterline. A synergistic retrieval algorithm has been developed to combine SMOS and CryoSat-2 data to estimate SIT more accurately (Kaleschke and others, Reference Kaleschke2015). This algorithm takes advantage of the strengths of both the data sets, incorporating SMOS thin ice information with CryoSat-2 freeboard measurements. By considering both thin and thick ice components, the algorithm can provide improved estimates of sea ice thickness across a wider range of ice conditions. The combined use of SMOS and CryoSat-2 data enhances the ability to monitor and study sea ice thickness variations over large spatial scales. By integrating the information from these two missions a more comprehensive view of the distribution and changes in SIT can be obtained. This large scale information is crucial for understanding the impacts of climate change over polar regions. This kind of SIT estimation algorithm is also used as part of the operational CMS service (Ricker and others, Reference Ricker2017) to provide weekly SIT over the Arctic Ocean. The SIT products are provided by the Sea Ice Thematic Assembly Centre (SITAC) operating under CMS.
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a remote sensing instrument that has also been utilized for thin sea ice thickness estimation (Yu and Rothrock, Reference Yu and Rothrock1996; Mäkynen and others, Reference Mäkynen, Cheng and Simila2013; Iwamoto and others, Reference Iwamoto, Ohshima and Tamura2014; Mäkynen and Karvonen, Reference Mäkynen and Karvonen2017). MODIS is aboard NASA's Terra and Aqua satellites and provides measurements in multiple spectral bands, including visible, near-infrared, and thermal infrared frequencies. While MODIS itself does not directly measure sea ice thickness, it offers valuable information for inferring sea ice characteristics. MODIS thermal infrared bands can provide measurements of sea ice surface temperature. Cold areas typically indicate the presence of thicker ice, as thick ice has a lower temperature compared to thinner ice. By examining the temperature gradient across the ice surface, ice thickness can be estimated. MODIS observations of surface reflectance can also be used to estimate ice surface albedo (Qu and others, Reference Qu2016). Albedo is related to SIT and the condition of sea ice (Barry and others, Reference Barry, Serreze, Maslanik and Preller1993), with thicker ice generally having a lower albedo due to its reduced light transmission. Thus, by analyzing changes in albedo, variations in sea ice thickness can be observed. It's important to note that while MODIS data can offer insights into sea ice thickness, they are often used in conjunction with other remote sensing techniques and field measurements to improve the SIT estimate accuracy and reliability. Also, use of other instruments operating at visible and near-infrared frequencies for sea ice thickness estimation has been studied, for example by Gu and others (Reference Gu2021). All instruments operating at visible and infrared frequencies have the problem that they are unable to operate in cloudy conditions. Utilization of visible frequency channels also requires daylight.
Synthetic Aperture Radar (SAR) is a widely used remote sensing technique intensively applied also in sea ice monitoring. SAR sensors, such as C-band SAR on RADARSAT-2, Sentinel-1 and Radarsat Constellation Mission (RCM), and L-band SAR on Japanese Space Agency's (JAXA) ALOS2/PALSAR2, emit microwave pulses and measure the backscattered signal from the sea ice surface and upper sea ice volume. SAR has proven to be valuable for estimating sea ice properties in a high resolution due to its all-weather and day-night independent imaging capabilities, thus enabling continuous monitoring of sea ice conditions. SAR resolution is typically from less than 10 m to a few hundred meters, depending on the instrument and its acquisition mode, i.e. significantly higher than typical passive microwave radiometer resolution ranging from a few kilometers to tens of kilometers. Existing X-band SAR instruments also have a good temporal and spatial coverage over the Baltic Sea because of their relatively wide swath width from about 100 km up to 270 km. The major areas of interest for sea ice monitoring in the Baltic Sea are the Gulf of Bothnia and the Gulf of Finland. These areas are typically covered by the current X-band SAR acquisitions within 1–2 days.
C-band SAR is more sensitive to sea ice with a thickness below 12 cm and L-band SAR frequency (1–2 GHz) penetrates deeper into sea ice and is sensitive to sea ice thickness between 12 cm and 36 cm. Because the X-band frequency (8–12 GHz) is higher than the C-band frequency (4–8 GHz) X-band SAR is sensitive only to even thinner sea ice thickness range than C-band SAR (Johansson and others, Reference Johansson, Brekke, Spreen and King2018). The study of the sensitivity of SAR instruments operating at different frequency bands was made for Arctic new ice. Because of lower salinity in the Baltic Sea we can expect similar or a little deeper penetration depth than observed in the above-mentioned study. Based on the penetration depth of typical SAR frequencies into sea ice, SAR data can only be used for thin ice thickness estimation directly. Because of the deeper penetration of L-band SAR it is the preferable SAR band for direct thin sea ice thickness estimation. C- and L-band SAR co-polarized channel backscattering ratio (VV/HH) for sea ice thickness estimation has been studied by Nakamura and others (Reference Nakamura2006) and the results indicate that backscattering ratios can be used to improve C- and L-band ice thickness estimation. The capability of C-band compact polarimetric SAR imagery has been studied in (Zhang and others, Reference Zhang, Dierking, Zhang, Meng and Lang2016). Some methods based on indirect inference with utilization of background SIT information have been proposed, e.g. in Karvonen and others (Reference Karvonen, Similä and Heiler2003, Reference Karvonen, Cheng and Similä2008). L-band SAR has been used for SIT estimation in Toyota and others (Reference Toyota, Ono, Cho and Ohshima2011). As a special case, estimation of pancake ice thickness based on SAR wave spectra has been studied in Wadhams and others (Reference Wadhams, Aulicino, Parmiggiani, Persson and Holt2018); De Carolis and others (Reference De Carolis, Olla and De Santi2021)
Also data fusion from SAR and other EO instruments acquired with a short time difference for SIT estimation has been studied. SIT estimation combining data from an ice model, MWR and SAR has been studied in Similä and others (Reference Similä, Mäkynen, Cheng and Rinne2013). Combining data from CryoSat-2 RA and SAR has been studied in Karvonen and others (Reference Karvonen, Rinne, Sallila, Uotila and Mäkynen2022). It seems that in the future the best SIT estimation results will be achieved by data fusion utilizing data from multiple EO instruments and ice models.
In this study we propose and evaluate a SIT estimation algorithm based on near-real-time X-band SAR data, received at FMI as part of CMS, and background SIT data derived from daily Baltic Sea ice charts and SIT grid provided by a thermodynamic ice model. This information is important for CMS and X-band SIT is needed to provide a better temporal SIT estimate coverage over the Baltic Sea. The temporal coverage of the C-band SAR SIT was significantly reduced after a critical anomaly of Sentinel-1B in December 2021. The anomaly ended the Sentinel-1B mission leaving only Sentinel-1A data available. In the current situation FMI is receiving a restricted amount of C-band data in a wide swath mode and HH/HV polarization combination from Radarsat-2 and Radarsat Constellation Mission (RCM). These C-band data are used for the operational SIT estimation but they do not provide an adequate temporal and spatial coverage over the Baltic Sea ice. For this reason, also X-band data are acquired over the Baltic Sea ice cover to complement the available C-band data. This study has been performed to be able to utilize the available X-band SAR data in SIT estimation and thus complement the temporally and spatially sparse C-band SIT currently provided for the CMS service.
2. Study area, data sets used in the study
The study area, Baltic Sea (see Fig. 1), is a semi-enclosed brackish water basin with a seasonal ice cover located in northern Europe, approximately between latitudes 53N and 65N and longitudes 9E and 30E. The area of Baltic Sea is 422 000 km2, and average annual maximum ice cover extent is 170 000 km2. Many of the harbors in the Baltic Sea are ice-surrounded every winter and precise and timely ice information is necessary for navigation in the ice-covered areas. The winter ship traffic is maintained with the aid of ice breakers. A typical Baltic Sea ice season lasts from November-December until late May in the northern parts (Gulf of Bothnia) and some weeks shorter in the Gulf of Finland. The thermodynamic ice growth in the fast ice zone is in maximum around one meter, on average 72 cm (Ronkainen, Reference Ronkainen2013), in the northern Gulf of Bothnia. In deformed ice areas ice thickness can be several meters, in ice ridges even 25 m (Kankaanpaa, Reference Kankaanpaa1997; Granskog and others, Reference Granskog, Kaartokallio, Kuosa, Thomas and Vainio2006). Maximum ice extent in Baltic Sea is typically reached in February-March.
Figure 1. Study area, locations of the SIT measurements (white dots) and areas of SAR image coverage (black frames).
The digitized FMI ice charts have been used in this study to provide daily SIT grid, i.e. minimum, average and maximum SIT at each grid point. In the ice charting procedure the ice parameters are estimated by ice analysts for polygons they draw on a map base according to their interpretation of the ice conditions. Each polygon represents an ice type or multiple ice types that can uniquely be described in terms of the ice charting guidelines provided by the World Meteorological Organization (WMO) (JCOMM, 2014; Canadian Ice Service, 2006). minimum, average and maximum SIT are also assigned to each polygon in the Baltic Sea ice charts. The FMI ice charts over the Baltic Sea are made by ice analysts daily during the winter period. The FMI ice charts are given in the Mercator projection. The same projection is used in the nautical charts over the Baltic Sea. The nominal resolution of the ice charts is about 1 km. The input data for making the FMI ice charts are satellite data from multiple instruments, including SAR (C band: Sentinel-1, RADARSAT-2, Radarsat Constellation Mission, and X band: COSMO-SkyMed, TerraSAR-X, PAZ), optical/infrared data (MODIS, VIIRS), observation data from coastal observers, observations from the Finnish and Swedish ice breakers and the FMI operational sea ice forecast model results. The most important data source are the Sentinel-1 C-band SAR images. In the FMI ice charts the ice analyst first locates the areas with homogeneous ice conditions. These areas are then presented by polygons in the ice charts. For each ice chart polygon five attributes (ice concentration, minimum ice thickness, average ice thickness, maximum ice thickness, degree of ice deformation) are assigned by the ice analyst. The ice concentration given in the FMI Baltic Sea ice charts is the total ice concentration of each ice chart polygon in percents. This differs from the WMO standard guidelines according to which sea ice concentration is given in tenths (in the included egg diagrams). It should also be noted that in the FMI Baltic Sea ice charts no partial ice concentrations neither stage of ice development are given. Instead, the five attributes listed above are assigned to each sea ice polygon. Over open water areas (polygons) the sea surface temperature (SST) attribute is also assigned instead of the sea ice attributes. The ice thickness values are level ice thicknesses, i.e. values for thermodynamically grown ice, possibly drifted from its origin. Ice deformation is given as a five-stage scale in which the value of one represents smooth level ice and the value of five highly deformed ice. The thickness of deformed ice can roughly be estimated by multiplying the level ice thickness by the degree of deformation. For example, rafted ice, with two overlying level ice layers, would then correspond to deformation degree of two.
The SAR data in this study were X-band data from COSMO-SkyMed (CSK), TerraSAR-X (TSX) and the Spanish PAZ X-band SAR (PAZ is peace in Spanish), all with HH polarization. CSK images were acquired in the ScanSAR Huge mode (200 × 200 km) and TSX and PAZ were acquired in the wide ScanSAR mode (270 × 200 km). PAZ SAR has similar properties as TerraSAR-X. During the season 2021–2022 400 X-band SAR images (216 TSX and 184 CSK images) were acquired over the Baltic Sea. In the training dataset we were able to use 100 TSX images and 91 CSK images of all the acquired images because only this fraction of the images could be associated with SIT measurements made during their acquisition date. During the season 2022–2023 totally 224 X-band SAR images (65 TSX images, 80 CSK images and 79 PAZ images) were acquired over the study area. However, spatially and temporally matching SIT measurements existed only for less than half of the images (totally 94 images: 32 TSX, 24 CSK and 38 PAZ images) and they were used in the test dataset.
The reference data were SIT measurements made by the Finnish and Swedish ice breakers along their routes in the ice covered Baltic Sea. The SIT measurements are point measurements. Because the ice thickness has some local variation and even large SIT differences within s small spatial range are possible the SIT measurements are typically rounded to the nearest 5 cm, depending on the measurer. Thus, we can assume that the accuracy of the ice breaker in-situ measurements is less than 5 cm. The SIT measurements performed during the SAR acquisition day were used in this study. Due to this selection the SIT measurement and the estimated SIT product do not represent exactly the same time instant. This may cause some additional estimation error with respect to the reference SIT data. Because of the very limited amount of measurements we had to make this compromise.
The FMI high-resolution thermodynamic snow and ice model (HIGHTSI) is used in this study to simulate thermodynamic ice growth and melt. HIGHTSI was initially developed to study snow and ice thermodynamics in seasonally ice-covered seas (Launiainen and Cheng, Reference Launiainen and Cheng1998; Vihma and others, Reference Vihma, Uotila, Cheng and Launiainen2002; Cheng and others, Reference Cheng, Launiainen and Vihma2003, Reference Cheng, Vihma, Pirazzini and Granskog2006). The model has been further developed to investigate snow and sea ice thermodynamics in the Arctic (Cheng and others Reference Cheng2008a, Reference Cheng, Vihma, Zhang, Li and Wu2008b). An operational version of HIGHTSI is run over the Baltic Sea daily during the Baltic Sea winter season.
Ice and snow thickness, heat conductivity and temperature are simulated, solving the heat conduction equations for multiple ice and snow layers. In addition, special attention is paid to the parametrization of the air–ice fluxes and the solar radiation penetrating into the snow and ice. The turbulent surface fluxes are parametrized taking the thermal stratification into account. The penetration of solar radiation into the snow and ice depends on the cloud cover, albedo, color of ice (blue or white), and optical properties of snow and ice. The shortwave radiation penetrating through the ice surface layer is parametrized to enable the model to calculate the sub-surface melting. Short- and long-wave radiative fluxes can either be parametrized or prescribed based on results of numerical weather prediction (NWP) models. The surface temperature is then solved from a detailed surface heat/mass balance equation, which is defined as the upper boundary condition and also used to determine whether surface melting occurs. At the lower boundary, the ice growth/melt is calculated on the basis of the difference between the ice–water heat flux and the conductive heat flux in the ice. The ice–water heat flux (refers to oceanic heat flux in this study) is assumed as a constant.
Model experiments in the Arctic Ocean have suggested that albedo schemes with sufficient complexity can more realistically reproduce basin-scale ice distributions (Liu and others, Reference Liu, Zhang, Inoue and Horton2007). We applied a sophisticated albedo parametrization according to temperature, snow and ice thickness, solar zenith angle and atmospheric properties (Briegleb and others, Reference Briegleb2004). The thermal conductivity of sea ice is parametrized according to Pringle and others (Reference Pringle, Eicken, Trodahl and Backstrom2007), and snow density and heat conductivity are calculated according to Semmler and others (Reference Semmler, Cheng, Yang and Rontu2012).
In seasonally ice-covered seas, such as Baltic Sea, snow plays an important role in sea ice growth. On one hand, snow generates a strong insulation and prevents sea ice growth. On the other hand, snow may transform to granular ice due to refreezing of ocean flooding and snow melt, enhancing ice growth. Snow-to-ice transformation through re-freezing of flooded or melted snow are considered in the model. A reliable calculation of snow depth heavily relies on precipitation, which is given as external forcing. HIGHTSI model over the whole Baltic Sea is run at FMI daily during the ice season and most recent available Baltic Sea HIGHTSI SIT grid before each X-band SAR acquisition, interpolated to the operational 500 m resolution from its rather coarse spatial resolution of approximately 15 km, has been used as one inputs in the proposed SIT estimation.
3. Methodology
The SAR images are first calibrated to logarithmic backscattering coefficients (σ 0) and after the calibration a linear incidence angle correction is applied to the calibrated backscattering coefficient (σ 0). The linear correction is performed such that σ 0 is corrected to correspond to the σ 0 of the incidence angle of 35 degrees. The slope used in the incidence angle correction is −0.25 dB/degree. This value was defined for the COSMO-SkyMED imagery acquired over the Baltic Sea during the winter season 2015–2016 to be able to correct the X-band SAR imagery incidence angle dependence in the operational SAR processing for the Finnish ice service. Totally 47 COSMO-SkyMED ScanSAR HUGE mode (frame size 200 km × 200 km) images were used in defining the incidence angle dependency. The dependency was defined by linear regression (a least squares fit). The slope was defined for sea ice areas only, distinguished from open water areas based on the corresponding day ice chart sea ice concentration. Open water areas were excluded because for open water the incidence angle dependence is dependent on the instantaneous local wave spectrum and over open water exact incidence angle correction without knowledge of the local instantaneous wave spectrum is impossible, and not even necessary for the scope of this study concentrating on sea ice only. The slope was rather similar as for the Baltic Sea C-band HH-polarized SAR imagery defined in Makynen and others (Reference Makynen, Manninen, Simila, Karvonen and Hallikainen2002). After SAR calibration and incidence angle correction, georectification to Mercator projection, also used in the Baltic Sea nautical charts, and land masking are performed, followed by the SAR segmentation. The land mask applied is based on the GSHHG (Global Self-consistent, Hierarchical, High-resolution Geography Database) coastline data (Wessel and Smith, Reference Wessel and Smith1996). We apply the iterated conditional modes (ICM) algorithm (Besag, Reference Besag1986) for the SAR segmentation in this study. ICM is currently used operationally at FMI for the CMS SAR-based sea ice products. SIT estimation is performed for each segment produced by the ICM segmentation, i.e. the SIT estimates are segment-wise. The features used in SIT estimation are median values within each segment. After preprocessing of the SAR image and the ICM segmentation, segmentwise features are extracted and they are then used in the segment-wise SIT estimation.
A general level diagram of the algorithm is presented in Fig. 2. After calibration and georectification some SAR texture features are computed and the ICM segmentation is performed. After this the segment-wise texture features and HIGHTSI SIT are computed as segment medians and based on these segment-wise values SIT estimates are computed. Finally, the segment-wise SIT estimates are scaled to be in the range of the FMI ice chart SIT based on the segment-wise minimum and maximum SIT from the digitized ice chart. The texture features were selected such that they have reasonably large absolute correlation value (>0.25) with SIT. We also tested some other texture measures, e.g local standard deviation but they were not included because of low their correlation with respect to SIT. The texture feature covariances were not taken into account in the selection of the texture features, and thus there exists at least some redundancy in the features.
Figure 2. General level diagram of the X-band SIT algorithm.
Because the amount of training data were very limited, consisting of totally 247 SIT measurements made by the Finnish and Swedish ice breakers during the winter 2021–2022, it was not possible to apply advanced machine learning (neural networks). Advanced machine learning algorithms require a large training data set. For this reason we here used a simple linear approach. The algorithm coefficients were defined by a linear regression based on least-squares fit applied to SAR backscattering coefficient (σ 0), four texture features derived from SAR imagery and the ice thickness provided by the FMI operational thermodynamic ice model HIGHTSI. The SIT estimate H at the 500 m grid location for each segment S, denoted by H S is computed in the first phase as:
where c i's and b j's are coefficients resulting from the linear regression and f i's are four texture features, explained in more detail later in this section, complemented by σ 0 and SIT values from the ice chart (average SIT) and HIGHTSI model. The second sum term in Eq. (1) represents the HIGHTSI model SIT modulated by σ 0 and the four texture features. The segment features used for a segment S are median values of the features within the particular segment.
100 CSK and 91 TSX images acquired during the winter 2021–2022 were used in the training. PAZ data for the season 2021–2022 were not available for training. Thus, we were not able to include PAZ in the training data.
The texture features, in addition to SAR σ 0, used were local autocorrelation, entropy, zero-crossing count and coefficient of variation (σ/μ).
Entropy E (Shannon, Reference Shannon1948) was computed as
where p k's are the proportions of each gray tone k within each sliding computation window. Auto-correlation C A (Box and Jenkins, Reference Box and Jenkins1976) 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 sliding computation window is here denoted by W. σ W and μ W are the mean and standard deviation within the window, respectively. The coefficient of variation is computed as
where σ and μ are the standard deviation and mean within the window w. The count of zero crossings is the number of separate connected areas with pixel values exceeding local pixel value average computed over each SAR segment normalized by the segment area:
where N ZC is the number of pixel blocks exceeding the segment average value withing a segment and A is the segment area (pixel count).
The significance of the features in the SIT estimation model can be estimated by the absolute value of the corresponding linear regression coefficient. Based on the regression coefficient magnitude the most significant features in the linear regression were ice chart SIT (H I), Z A, HIGHTSI SIT (H H) and σ 0 but also the other included texture features had some significance. Of the modulation terms the product of H H and Z A was the most significant term. We do not provide the regression coefficient values here because they will be frequently updated as more training data will appear during the Baltic Sea winter seasons. Redefining the coefficients is a simple least-squares fit and it is fast to perform to include the moderate amount of new SIT measurements e.g. annually after each winter season.
Because applying the proposed linear regression directly seemed to produce too thick sea ice, even twice as thick as the measured SIT at some locations, a restricting normalization was applied to the SIT produced by the regression as a post-processing phase. The normalization maps the thickness within a mapped ice chart segment to be between the minimum and maximum ice chart SIT within the mapped ice chart segment. The mapping is linear, such that the lowest estimated SIT is mapped to the minimum ice chart SIT and highest estimated SIT to the maximum ice chart SIT.
4. Results
Comparison between the estimated SIT and reference SIT measurement performed by the ice breakers during the season 2022–2023 were made by using four different measures, L1 difference (L1D) the signed L1 difference (bias, B) root mean square (RMS) difference RMSD and correlation (R):
N s is the number of evaluation SIT measurements $H_i^{\hbox{est}}$
is the estimated SIT, and $H_i^{\hbox{ref}}$
is the measured reference SIT at each pixel indicated by the index i. $\overline {H}^{\hbox{est}}$
and $\overline {H}^{\hbox{ref}}$
are the means of the estimated and reference SIT, respectively. L1D ( ≥ 0) describes the mean absolute error. It is a generally used error measure providing a good overview of the estimation accuracy. B describes the bias, i.e. the systematic error. If B < 0, the algorithm underestimates, and if B > 0, the algorithm overestimates the SIT with respect to the reference data. RMSD was included because larger differences have large effect on RMSD and thus RMSD is sensitive to outliers. Correlation R is a statistical measure of how the two variables, in this case H est and H ref, vary together, i.e. how well it is possible to predict the value of the other variable when the other variable value is known.
In Table 1 also the corresponding statistics for the ice chart (IC) SIT, HIGHTSI model SIT, SIT estimated by a linear combination of ice charts SIT and HIGHTSI SIT (IC+HIGHTSI), and operational C-band SIT estimation algorithm using Sentinel-1 (S-1) and Radarsat-2 (RS-2) data are given. The same reference data points were used in all the comparisons. The C-band SAR SIT used as an additional reference dataset in the table is produced by the FMI operational CMS SIT product (Karvonen and others, Reference Karvonen, Similä and Heiler2003). It can also be seen that the error statistics are rather similar for all the three studied X-band SAR instruments. Bias for them was very small, in maximum 1.5 cm, L1 difference on average around 7 cm and RMSD some over 9 cm. Ice chart SIT and C-band SAR SIT have small negative bias values (underestimation) but HIGHTSI has a large positive bias, i.e. HIGHTSI on average seems to overestimate the ice thickness at the SIT measurement points. The large bias can at least partly be explained by the fact that HIGHTSI is a thermodynamic ice model and it does not take the ice dynamics into account. For this reason HIGHTSI tends to overestimate SIT in areas where ice has drifted away and redeveloped causing the positive bias. RMSD is also larger for C-band and ice chart SIT, indicating that there are more large differences (outliers) compared to X-band SIT. For HIGHTSI SIT part of the large RMSD can be explained by the large bias. However, because of the coarse resolution and exclusion of ice dynamics in the model also large local differences exist. The correlations were moderate, over or very close to 0.5, for all the other SIT datasets except for HIGHTSI (0.3). SIT estimated as a linear combination of ice chart SIT and HIGHTSI SIT has been included in the comparison to exclude all the SAR data information from the SIT estimate. The linear model coefficients were defined based on the season 2021–2022 (training) data. The X-band SAR algorithm yielded significantly better results than the linear model excluding all SAR data. However, the performance measures of the linear model using ice chart and HIGHTSI SIT as its input were quite close to the results of the operational C-band SAR SIT algorithm. It should be noted that in the C-band SIT algorithm ice model data are not included.
Table 1. Difference measures between the ice breaker measurements and X-band SIT estimates for the season 2022–2023
Positive bias indicates overestimation. For reference, ice chart (IC) SIT, HIGHTSI SIT, SIT estimated as a linear combination of ice chart and HIGHTSI SIT and the operational Sentinel-1/Radarsat-2 C-band SAR SIT at the same reference points have been included.
In Fig. 3 the estimated X-band SIT is plotted against the ice breaker SIT measurements. It can be seen that in general. there exists some deviation, but because the in-situ measurements are point measurements and the SAR SIT values always present an average over a larger areas, in this case over segments of variable sizes, there will always exist some difference between SAR SIT and reference SIT measurements because the in-situ and SAR based SIT values do not describe exactly the same thing.
Figure 3. X-band SIT vs all the 2023 ice breaker measurements. The corresponding Pearson correlation coefficient value was 0.61.
In Figs. 4 and 5 examples of the X-band SIT estimation are shown. In Fig. 4 the HIGHTSI thermodynamic ice model SIT, the ice chart SIT, combined (linear model) ice charts and HIGHTSI SIT, an X-band SAR mosaic of 1 March 2023 and the corresponding X-band SAR SIT estimation are shown. The SAR imagery covers only the ice covered areas of the Baltic Sea, i.e. the Gulf of Bothnia and Gulf of Finland, the dark area in the SAR mosaic of Fig. 4d is the area outside of the acquired SAR imagery. In Fig. 5 HIGHTSI ice model SIT, ice chart SIT, combined ice chart and HIGHTSI SIT, and X-band SIT corresponding to a X-band SAR frame of 23 January 2024 over the Bay of Bothnia (northern part of Gulf of Bothnia) is shown. This X-bad SAR SIT example is taken directly from the operational CMS processing chain. It can be seen that the ice model SIT does not show much variation over the area. This is because the model is thermodynamic and the winter weather in December 2023 and January 2024 was rather cold and similar over the whole image area. The X-band SAR SIT exhibits much more details than the ice model SIT, the ice chart SIT or the combined ice chart and ice model SIT. These details are very useful for ice breaking and ice navigation enabling optimal route selection in sea ice.
Figure 4. X-band SIT estimation over the Baltic Sea on 1 March 2023: HIGHTSI segment median (a), ice chart SIT (b), linear model combining ice chart and HIGHTSI SIT (c), X-band SAR mosaic (d), and X-band SIT mosaic (e).
Figure 5. An example of X-band SIT over the Gulf of Bothnia, based on a single SAR image, 23 January 2024: HIGHTSI SIT segment-wise median (a), ice chart SIT (b), linear model combining ice chart and HIGHTSI SIT (c), and X-band SAR SIT (d).
5. Discussion and conclusions
We have introduced a SIT estimation algorithm based on HH-polarized X-band SAR data and background SIT information from the FMI ice charts and thermodynamic ice model here. The algorithm is applied as part of the operational CMS service and the operational X-band SIT service was started in December 2023. The CMS sea ice products are provided in NRT. The requirement for the SAR data is to be available from a ground station in less than three hours after the SAR acquisition and the CMS sea ice products to be available within four hours after the SAR acquisition, i.e. leaving one hour for data transmission from the ground station to FMI, algorithm execution and uploading the product to the CMS service. The most visible difference to the ice chart SIT and the modeled SIT is the higher level of detail in the SAR-based SIT. Because the SAR segments are more detailed and also their boundaries are more precise when compared to hand-drawn ice chart polygons, it is possible to distinguish the areas of different ice type and thickness with a higher precision and level of detail than available in the ice charts by combining SAR, ice model and ice chart information as proposed here. Inclusion of X-band SIT in the CMS product portfolio is important because of the reduced wide swath C-band SAR data in HH-polarization after losing Sentinel-1B. X-band SIT is able to cover the temporal and spatial gaps in the C-band SIT. From the beginning of the winter season 2024–2025 FMI wil also be providing a daily SIT mosaic product over the whole Baltic Sea combining the single C-band and X-band products for CMS.
The selection of texture features is now based on the correlation magnitude between the measured SIT and the features. This does not take into account the inter-feature dependencies (covariances). For example, the Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, Reference Tibshirani1996) could be used to select the features among a set of texture features taking the inter-feature dependencies into account. However, the current set of texture features seems to work well in operation and the next generation algorithm will very likely utilize a more sophisticated machine learning approach instead of a simple linear regression, and only minor updates to the current algorithm will be made before moving to a next generation algorithm. We also included so-called modulation terms combining the HIGHTSI SIT and the SAR features by multiplication (HIGHTSI SIT modulated by the SAR features) in the algorithm. These terms were included to involve local variation in the low-resolution HIGHTSI SIT, alone representing average coarse-scale thermodynamic SIT without ice dynamics taken into account. Linear regression was used here because of the small amount of available SIT measurements. We are studying possibilities to synthesize large SIT training datasets from the original SIT measurements and auxiliary information from model and EO data to be able to train more advanced machine learning algorithms for SIT estimation in the future.
According to the evaluation against the ice breaker SIT measurements the accuracy of the X-band SIT product was at an acceptable level for the operational use. The evaluation statistics using the same evaluation data were better than for the current FMI operational CMS SITAC SIT products, produced based on C-band SAR data from Sentinel-1 and Radarsat-2. We at FMI have started development of a new C-band SIT algorithm utilizing both HH and HV polarization channels available in the C-band imagery. The new C-band algorithm will utilize texture features computed for both the polarization channels and the HH/HV channel cross-correlation. The development of this new C-band SIT algorithm is in a very preliminary phase and reliable performance measures can be provided in a later phase of the development. The X-band SIT RMSD values were relatively high compared to the L1D values. This indicated that there were some larger single differences with respect to the reference SIT measurements. One reason for these large differences is the different nature of the measurements and estimates, i.e. the measurements are point measurements and the estimates represent averages of larger areas (segments). The ice charts SIT is given for even larger areas, described by the ice chart polygons. This difference in the nature of the SIT measurements and estimates explains at least part of the larger RMSD of the ice chart SIT and the HIGHTSI SIT when compared to the X-band SAR SIT. The ice thickness may vary locally within a short distance (even within only a few meters). Because of the small amount of SIT measurements it was impossible to estimate e.g. segment SIT averages because there were only single measurement points for a single segments both in the training and test datasets. The SIT measurements were measured during the SAR acquisition day and thus there typically exists a time gap of several hours between the SAR acquisition and the SIT measurement. This will increase the estimation error at least in some cases when the ice has been moving rapidly. However, as the SIT estimates were provided for segments representing a larger area than single points the estimation is more robust to moderate ice motion between the measurements time and SAR acquisition time than point-wise estimation. The reference measurements are mostly made in the drift ice field and it seems that the algorithm performs well on both landfast ice zone and drift ice zone. Actually, in the landfast ice zone the average ice thickness can rather well be estimated based on the thermodynamics alone, e.g. by the HIGHTSI model. In the drift ice zone the role of SAR imagery becomes more important. SAR backscattering and texture indicates the locations of thin ice (typically smooth level ice) and thicker deformed ice (rough surface). This is essential information provided by SAR for the SIT estimation.
The algorithm was run in an operational test mode during the season 2022–2023 and is now, starting from the beginning of December 2023, run in the fully operational mode as part of the CMS SITAC Baltic Sea ice product portfolio. The CMS products will be evaluated after each season in May-June and the seasonal evaluation against the SIT measurements will be used for tuning the algorithm after each season. After running the algorithm for a few seasons more reference SIT measurements will be available. The increased amount of reference data can be used to improve the algorithm, for example by testing the use of different parametrization for different phases of winter. For example, early winter freeze-up, steady mid-winter and late winter/spring meltdown phases could have different linear coefficients produced by separate linear regression fits.
One rather straightforward future evolution step will be to derive the range of sea ice thickness from the modeled sea ice (HIGHTSI) and completely exclude the manually produced ice chart information from the processing. The local SIT minimum and maximum could, for example, be the minimum and maximum of the modeled SIT within a given range from each SIT estimation grid cell location. This would enable a SIT product that is independent of the manual ice charts. Also, methods to derive large synthetic SIT datasets for training of more sophisticated machine learning based SIT estimation models will be studied. Generation of the synthetic data sets could utilize the available in-situ measurements and additionally all the available SIT related information from EO data. This would be an important step towards an advanced machine learning SIT estimation algorithm involving SAR data, model data and other available EO data, e.g. from altimeters and microwave radiometers.
Open access
