An Arctic sea ice concentration data record on a 6.25 km polar stereographic grid from 3 years of Landsat-8 imagery

The decline in Arctic sea ice in the global warming era has received much attention as a contributing factor to the changes in the weather and climate in the Arctic and beyond. The coverage of Arctic sea ice (i.e. sea ice concentration (SIC)) has been monitored since 1972 using satellite passive microwave (PMW) measurements because of their extensive coverage and all-weather capability. However, the fundamental basis of algorithms for estimating SIC has not improved much since the early days due to the lack of reference SIC data, leading to discrepancies between existing PMW SIC algorithms. To overcome this issue, this study aims to construct data records of reference SIC over Arctic sea ice using 30 m resolution imagery from the Operational Land Imager (OLI) on board Landsat-8. In order to collect relatively bright and clear scenes, thresholds of solar elevation >15° and cloud cover <10 % were applied in this study. Clouds in each Landsat-8 scene were masked using the cloud-masking array provided in Landsat data. Due to the poor accuracy of the cloud-masking array over ice-covered surface types, an additional step of visually inspecting the state of the cloud mask using the true-colour image was designated in this study. Each Landsat-8 scene was sorted into four categories depending on the state of the cloud mask. The normalized difference snow index and OLI band-5 reflectivity were used to differentiate between ice and open water within each selected Landsat-8 pixel. The classified data were projected onto a 6.25 km polar stereographic grid, and SIC for each grid cell was obtained by counting ice-classified pixels within the grid. SIC was only computed for grid cells where more than 99 % of their area was covered with Landsat-8 pixels to limit the uncertainty in SIC arising from grids that are not fully concentrated with Landsat-8 pixels. Uncertainty in the produced SIC was 1 %–4 %, inferred using the Gaussian error propagation method. Out of 15 286 collected Landsat-8 images, 14 297 images were translated into SIC maps, and a total of 2 934 399 Landsat-8 SIC grid cells were generated. Evaluation of Landsat-8 SIC with SIC from ice charts revealed a good linear relationship (correlation coefficient of 0.96) between the two products, as well as a mean negative bias which fell within the uncertainty range of Landsat-8 SIC. SIC based on Landsat-8 can be used as reference SIC to evaluate existing SIC products, and, thus, one can improve SIC products, as well as use the improved SIC for other applications such as data assimilation and retrieval studies. The vast amount of Landsat-8 SIC generated in this study may also be used to train deep-learning models for the estimation of Arctic SIC coverage. The Landsat-8 SIC dataset can be publicly accessed at https://doi.org/10.5281/zenodo.10973297 (Jung et al., 2024), and the Python codes for the production and evaluation of the Landsat-8 SIC dataset are accessible at https://doi.org/10.5281/zenodo.12754602 (Jung, 2024).

1 Introduction

Since space-borne multi-channel passive microwave (PMW) observations have been available, areal information of Arctic sea ice has been successfully monitored. During the past 4 decades, these observations have shown that sea ice extent (SIE), which is defined as the area of ocean where sea ice concentration (SIC) is greater than 15 %, has been rapidly declining at a statistically significant negative trend of −12.7 % per decade, as observed in September (Cavalieri and Parkinson, 2012; Meier et al., 2014; Meier and Stroeve, 2022). In the global warming era, the change in Arctic sea ice area is considered to be a key indicator of climate change, and this is closely associated with the changes in the Arctic local weather, as well as in the weather at mid-latitudes (Honda et al., 2009; Jaiser et al., 2012; Kim et al., 2014; Trewin et al., 2021; Shi et al., 2023). Therefore, obtaining precise observations of Arctic SIC is essential in order to identify the influences of climate change on Arctic sea ice.

As mentioned above, the spatial coverage of Arctic sea ice (i.e. SIC) has been monitored using satellite PMW measurements, with their extensive spatial coverage over the Arctic and all-weather capability. Beginning with the launch of the Electrically Scanning Microwave Radiometer (ESMR) on board Nimbus-5 in 1972 (Parkinson et al., 1987), successive launches of PMW sensors have allowed for the construction of comprehensive and continuous records of Arctic SIC. The Scanning Multichannel Microwave Radiometer (SMMR), launched in 1978, was equipped with five channels (6.6, 10.7, 18.0, 21.0, and 37.0 GHz) in the first two Stokes' polarizations. The emergence of multi-channel PMW radiometers has led to the development of various SIC retrieval methods which were more accurate relative to the previous methods used for the ESMR, which only had a single channel at 19 GHz. The addition of the near-90 GHz high-frequency channels into the PMW sensors following the SMMR, which include the Special Sensor Microwave Imager (SSM/I), the Special Sensor Microwave Imager/Sounder (SSMIS), the Advanced Microwave Scanning Radiometer – Earth Observing System (AMSR-E), and the AMSR2, has allowed for spatially enhanced SIC retrievals.

Various PMW SIC algorithms have been developed which estimate SIC based on combinations of brightness temperatures (TBs) at various channels and empirically derived tie points. One of the best-known algorithms is the bootstrap (BT) algorithm first suggested by Comiso et al. (1984). In the BT algorithm, vertically polarized TBs at 19 and 37 GHz and horizontally polarized TB at 37 GHz are utilized to determine reference TBs (i.e. tie points) over open water and fully concentrated ice, which can be used to convert the observed TB to SIC with the following equation:

$$\begin{array}{}\text{(1)}& \mathrm{SIC}={\displaystyle \frac{T_{\mathrm{B}}-T_{\mathrm{O}}}{T_{\mathrm{I}}-T_{\mathrm{O}}}},\end{array}$$

where T_B is the satellite-measured TB, and T_O and T_I are the empirically determined open-water and ice tie points, respectively. Tie points in the BT algorithm are updated on a daily basis and are acquired separately for the Arctic and the Antarctic in order to accommodate the variation in TB fields with respect to time and hemisphere (Comiso, 1995). Another well-known algorithm is the NASA Team (NT) algorithm, which utilizes horizontally polarized TB at 19 GHz and vertically polarized TBs at 19 and 37 GHz to calculate the polarization ratio (PR) and the spectral gradient ratio (GR), which are used to determine a set of tie points to estimate SIC and to determine the surface type from a combination of open water, first-year ice, and multi-year ice (Cavalieri et al., 1984). The other is the Arctic Radiation and Turbulence Interaction Study (ARTIST) Sea Ice (ASI) algorithm, which was developed by Kaleschke et al. (2001) in order to exploit the high resolution of the near-90 GHz channels. The ASI algorithm estimates SIC using the tie-points derived from the polarization difference calculated in the near-90 GHz channels. The high sensitivity of the near-90 GHz channels to atmospheric effects is compensated for by the usage of weather filters, which are applied using the GR thresholds suggested by Gloersen and Cavalieri (1986) and Cavalieri et al. (1995) and by setting SIC to zero in areas where BT SIC values are zero (Spreen et al., 2008).

However, discrepancies exist among various PMW SIC records retrieved from different algorithms owing to the different channel combinations, tie points, and weather filters used in each algorithm (Comiso et al., 1997; Andersen et al., 2007). Due to the lack of reference SIC data with satisfactory temporal and spatial coverages, these disagreements have been studied mainly through the intercomparison of different PMW SIC values and ensemble methods which compare individual SIC products to their averaged value. For instance, Ivanova et al. (2014) reported that different PMW SIC products showed a maximum difference of up to 1.3 × 10⁶ km² in area and 0.6 × 10⁶ km² in extent over the Arctic and larger deviations during the summer due to the differing sensitivity of retrieval algorithms to the presence of melt ponds and the associated emissivity change, as well as to a humid atmosphere (Ivanova et al., 2014; Comiso et al., 2017; Horvat et al., 2023). Although these intercomparison approaches can provide valuable assessments of the consistency of PMW SIC products from sub-seasonal to climatological timescales, it is noted that there is a limitation in providing a quantitative assessment of PMW SIC products.

In order to make such quantitative assessments, it is essential to have independent SIC data that can be used as a reference. Spaceborne sensors with visible (VIS) to infrared (IR) channels, such as the Moderate Resolution Imaging Spectroradiometer (MODIS) and the sensors on board the Landsat series, have been used to generate reference SIC due to their finer spatial resolutions compared to PMW sensors (Markus et al., 2002, Cavalieri et al., 2006, 2010; Rösel and Kaleschke, 2011; Kern et al., 2022; Tanaka and Lu, 2023; Song and Minnett, 2024). However, validation attempts using VIS- and/or IR-based SIC as a reference have been limited to the usage of a small number of VIS and/or IR images, with the exception of the dataset of Kern et al. (2022), which used a relatively large number of Landsat scenes (386 scenes) to generate a reference SIC. The dataset by Kern et al. (2022) is also utilized in this study for validation of the produced Landsat-8 SIC, and the results of the comparison are presented in Sect. 3.2. In addition to VIS and/or IR instruments, SIC observations from synthetic aperture radar (SAR) have also been used for PMW-based SIC validation purposes, but difficulties in obtaining an accurate and automated SIC map from SAR images result in the limited use of SAR images for validation purposes (Andersen et al., 2007; Han and Kim, 2018; Park et al., 2020; Tanaka and Lu, 2023).

Recently, efforts to leverage the advantages of both VIS and/or IR sensors and PMW sensors for retrieving SIC have been explored through data-merging techniques. Ludwig et al. (2020) used a combination of MODIS and AMSR2 measurements to construct a high-resolution (1 km) and spatially continuous SIC dataset over pan-Arctic areas. This approach exploited the benefits of the 1 km resolution MODIS imagery while mitigating its inherent disadvantage of spatial discontinuity due to clouds by introducing the AMSR2 measurements. While the SIC dataset produced by Ludwig et al. (2020) is both high-resolution and covers pan-Arctic areas, due to the retrievals being reliant on the AMSR2 measurements, the product cannot be considered to be a fully independent reference dataset for PMW SIC validation. Therefore, it is still necessary to construct a dataset of Arctic SIC that is fully independent of PMW measurements.

In addition to this, recent applications of deep-learning (DL) models for estimating SIC have shown promising results. Karvonen (2017) trained a multi-layer perceptron (MLP) model using various combinations of PMW signals extracted from AMSR2 and SAR as the training inputs and SIC fields derived from the Finish Meteorological Institute ice charts as the reference. This MLP model produced improved SIC compared to the high-resolution ASI SIC. However, the data used to train the DL model suggested by Karvonen (2017) were limited to regions around the Baltic Sea and were only acquired during the winter of 2015–2016. Chi et al. (2019) proposed an estimation of Arctic SIC based on an MLP model trained with raw AMSR2 TBs as the inputs and SIC derived from 72 MODIS images during 2016 as the reference, demonstrating that the DL-based SIC shows better performance than the widely used BT and ASI SIC. Since both studies used training datasets acquired during a limited time period but showed promising results in terms of the use of DL techniques for SIC production, it is also desirable to construct a data record for reference SIC data with satisfactory temporal and spatial coverages.

Therefore, this study aims to construct a reference SIC dataset of satisfactory spatiotemporal extent that allows for the validation of various SIC products over pan-Arctic areas and that can be used for DL training. To do this, a total of 14 297 Landsat-8 images over 3 years (2020–2022) were translated into SIC maps in a 6.25 km polar stereographic grid and were catalogued into a region of the Arctic Ocean.

The remaining sections of this paper are organized as follows: Sect. 2 provides a detailed description of the Landsat-8 dataset, the land, the sea ice region, the coastal area masks, and the reference datasets used to evaluate the Landsat-8 SIC in this study. Section 3 describes the pipeline of processing a Landsat-8 image into a SIC dataset, along with the uncertainty estimation. The resultant SIC product and its uncertainty are shown in Sect. 4. Possible errors in Landsat-8 SIC, evaluation of Landsat-8 SIC using existing SIC from ice charts, evaluation of Landsat-8 SIC over melt ponds, and qualitative assessment of two PMW SIC products using SIC from Landsat-8 as a reference are discussed in Sect. 5. Section 6 provides the data availability statement, and Sect. 7 presents the summary and conclusions of this research.

2 Used data

2.1 Landsat-8 OLI-TIRS Collection 2 Level 1 products

In this study, reflectivities measured by the Operational Land Imager (OLI) on board Landsat-8, which is a polar-orbiting satellite with an orbit inclination of 98.2° and a repeat cycle of 16 d (Zanter, 2019), were used to retrieve SIC values over pan-Arctic areas. The OLI sensor has a swath width of 185 km, measuring radiances in eight bands from VIS to shortwave IR (SWIR), with a spatial resolution of 30 m. It should be noted that areas with a latitude higher than 82° N in the Northern Hemisphere are not measured by Landsat-8 (i.e. the hatched area in Fig. 1) due to the orbit inclination of Landsat-8 and the relatively narrow swath width of the OLI. The Landsat-8 Collection 2 Level 1 product used in this study contains 11 spectral-band images (nine bands from the OLI and two bands from the Thermal Infrared Sensor (TIRS), provided in GeoTIFF format, and two quality assessment bands containing masking information for clouds, cloud shadows, cirrus, fill values, and radiometric saturation). To calculate the SIC, the OLI-measured reflectivities at near-infrared (NIR) band 5 and SWIR band 6 (used in the normalized difference snow index) were used in this study. It is worth noting that the methods developed in this study (described in Sect. 3) utilize the NIR and SWIR bands for SIC retrieval and are therefore applicable to a wider range of high-resolution sensors that observe at similar bands, including the Multi-Spectral Instrument (MSI) on board Sentinel-2. However, due to the more robust cloud-masking performance of the Landsat-8 product, in this study, the Landsat-8 Collection 2 Level 1 product was selected to be used for the production of reference SIC data (Zhu et al., 2015; Tarrio et al., 2020).

For the period of January 2020–December 2022, the Landsat-8 Collection 2 Level 1 product and the corresponding true-colour images were downloaded from the United States Geological Survey Earth Explorer (https://earthexplorer.usgs.gov/, last access: 16 March 2024). To circumvent the influences of cloud contamination and low solar elevation angle on SIC calculation, which hamper accurate classification of ice and open water, only Landsat-8 images with less than 10 % cloud cover (based on fractional cloud-masked area from the quality assessment band of Landsat-8) during the daytime (solar elevation higher than 15°) were collected. While the threshold of 0 % cloud cover would ensure the acquisition of the least cloudy scenes, this also results in the loss of a considerable number of Landsat-8 scenes that contain clear-sky portions (see Fig. S1 and Table S1 in the Supplement for the number of available Landsat-8 scenes subject to different threshold values of cloud cover). Therefore, the threshold value for cloud cover was relaxed to 10 % during the acquisition of Landsat-8 images. Since VIS measurements are not available during polar nighttime, the Landsat-8 data between early December to January were not collected. A total of 15 286 images were collected and sorted into 12 regions of the pan-Arctic area for the calculation of SIC. In the case of a Landsat-8 image that was observed across more than one region, the image was sampled repeatedly for each region. Footprints of the collected Landsat-8 images are displayed in Fig. 1, and the number and the temporal availability of the collected images for each area are listed in Table 1.

Table 1 The number of Landsat-8 images collected in this study and the available period for each region of the pan-Arctic areas.

Figure 1 Footprints of the collected Landsat-8 images over each region of the pan-Arctic areas during the period of January 2020–December 2022. The hatched region denotes the areas unmeasured by Landsat-8 due to its orbital inclination (i.e. pole hole). The regions of the pan-Arctic areas were distinguished using the region mask provided by Meier and Stewart (2023a). The map projection is the NSIDC Sea Ice Polar Stereographic North (EPSG: 3413), and the map was plotted using Python.

2.2 Land and sea ice region masks

Regions of the Arctic Ocean were distinguished using the National Snow and Ice Data Center (NSIDC) “Arctic and Antarctic Regional Masks for Sea Ice and Related Data Products, Version 1” data (Meier and Stewart, 2023a), which divide the Arctic Ocean into 19 different regions with 6.25, 12.5, and 25 km resolution polar stereographic (PSR) grids. In addition, this product provides surface-masking information to differentiate between ocean areas and non-ocean areas such as land, freshwater, land ice, ice shelf, and disconnected ocean. The flag values for the Arctic Ocean regions and the different surface types can be found in the product user guide (Meier and Stewart, 2023b). Amongst the regions, 12 regions (i.e. Baffin Bay, Barents Sea, Beaufort Sea, Bering Sea, Canadian Archipelago, Central Arctic, Chukchi Sea, East Greenland Sea, East Siberian Sea, Hudson Bay, Kara Sea, and Laptev Sea; Fig. 2) were selected to generate Landsat-8-based SIC values because the above 12 regions have climatologically meaningful sea ice extent.

Figure 2 Geographic distribution of the designated regions of the Arctic Ocean based on NSIDC Sea Ice Region Mask data (Meier and Stewart, 2023a). The map projection is NSIDC Sea Ice Polar Stereographic North (EPSG: 3413), and the map was plotted using Python.

2.3 Ice- or water-classified Landsat-8 images

In this study, the performance of ice and open-water classification (later described in Sect. 3.2) was evaluated using ice–water classification data from the “Land surface type over water from supervised classification of surface broadband albedo estimates” (Kern, 2021; Kern et al., 2022). This dataset contains ice and water classification estimates using broadband albedo values from the Landsat series (i.e. bandwidth-weighted mean albedo from Landsat-8-measured reflectivities at bands 3, 4, and 5), where each pixel in a scene is classified into open water, thin or bare ice, and thick or snow-covered ice based on supervised classification. In our study, the two ice categories (i.e. one for thin or bare ice and the other for thick or snow-covered ice) were considered to be the same ice category due to the higher ambiguities in the discrimination among different ice types relative to the discrimination between ice and open water (Kern et al., 2022). In order to evaluate the classification method suggested by our study, we processed Landsat-8 reflectance from six clear-sky scenes that Kern (2021) had classified and then compared results. The result of the comparison is presented in Sect. 3.2, and the location and time of the Landsat-8 scenes that were used in the evaluation are provided in Fig. S2 and Table S2 in the Supplement.

2.4 Ice chart data

Ice charts provide SIC intervals over the Arctic obtained by means of manual interpretation of satellite images from various sensors such as SAR, MODIS, and the Advanced Very High Resolution Radiometer (AVHRR). In this study, operational ice charts from the Norwegian Meteorological Institute (MET Norway), which provide SIC maps in a PSR grid with a nominal resolution of 1 km, were used to evaluate the performance of the produced Landsat-8 SIC. Each grid in the ice chart contains the six classified SIC values (5 %, 20 %, 50 %, 75 %, 95 %, and 100 %), which represent the ice concentration intervals defined by the World Meteorological Organization (WMO) (Table A1 in the Appendix). The ice charts are provided on a daily basis and cover the spatial domain of approximately 60–85° N, 80° W–80° E, which overlaps with the regions of the Barents Sea, Central Arctic, East Greenland Sea, and Kara Sea defined in Sect. 2.2. It is noted that SIC values in ice charts are based on the interpretation of multiple satellite images by ice analysts and therefore contain high uncertainties, which are reflected by the wide ice concentration intervals designated for each of the six SIC values (Table A1). Even with such high uncertainties, SIC values from ice charts have been widely selected as reference data in SIC product validation studies because they can be used to provide quantitative information about the observed ice coverage (Agnew and Howell, 2010; Ivanova et al., 2015; Karvonen, 2017).

In this study, 2 years (2021 and 2022) of ice charts were collected; among these, 222 ice charts that presented a spatial overlap with the coverage of Landsat-8 SIC and a time difference of less than 1 h compared to the Landsat-8 scene were used for the evaluation of the produced Landsat-8 SIC (see Table S3 in the Supplement for the list of ice charts and Landsat-8 files used in the evaluation of the produced Landsat-8 SIC).

2.5 Melt pond fraction data

Melt ponds are formed from the surface melting of sea ice and are known to exist in preponderance over the Arctic during the melting season (Untersteiner, 1961; Fetterer and Untersteiner, 1998; Rösel et al., 2012). In the VIS and IR ranges, melt ponds typically exhibit lower spectral reflectivities relative to dry sea ice (Perovich, 1996; Malinka et al., 2018) and therefore may introduce errors into SIC estimated from VIS and/or IR observations because the optical characteristics of melt ponds may not be differentiated from those of open ocean. In order to test the sensitivity of Landsat-8 SIC values to the existence of melt ponds, in this study, a melt pond fraction (i.e. the fractional areal coverage of melt ponds over sea ice; MPF) dataset estimated from clear-sky Sentinel-2 satellite imagery was introduced (Niehaus and Spreen, 2022; Niehaus et al., 2023). This dataset also contains an open-water mask (OW mask), which is a binary classification mask of each pixel in a Sentinel-2 scene into ice and open water. This dataset is available from 2017 to 2021 for the Arctic melting season (i.e. June, July, and August). In this study, each MPF dataset was tested for spatiotemporal overlap (time difference of less than 3 h) with the coverage of Landsat-8 SIC. A total of six MPF datasets were found to overlap with the coverage of Landsat-8 SIC and were thus available for use in the evaluation. The list of available MPF datasets and the corresponding Landsat-8 scenes can be seen in Table S4 of the Supplement.

3 Method

Figure 3 shows the pipeline of processing a level-1 Landsat-8 image into an SIC product based on 6.25 km resolution PSR grid. The details of each process are explained in the following subsections.

Figure 3 Pipeline for the processing of level-1 Landsat-8 OLI images into SIC with 6.25 km resolution.

3.1 Removal of cloud pixels and cloud mask quality assessment

Satellite observations of surface properties from the VIS and NIR channels are hindered by the presence of clouds. Therefore, it is important to filter the presence of clouds prior to the SIC data production. In this study, clouds and cloud shadows within each Landsat-8 scene were masked using the masking array constructed from the quality assessment band of each Landsat-8 level-1 product, which is generated by the C Function of Mask (CFMask) (Zhu and Woodcock, 2012). CFMask is a cloud detection algorithm that provides masking information for clouds, cloud shadows, and cirrus. Confidence scores are also given in three levels (i.e. low, medium, and high) for clouds and two levels (i.e. low and high) for cirrus. A confidence score for cloud shadows is not provided because cloud shadows are only derived from high-confidence cloud pixels by using the geometric relationship between the position of the sun and high-confidence cloud pixels (Zhu and Woodcock, 2012). Although the application of the lowest confidence scores in the removal of clouds and cirrus would ensure the lowest rate of false negatives (FNs; cloud pixels that are mistaken as clear pixels) in cloud detection, the use of the lowest confidence scores would also result in the removal of a considerable number of sea ice pixels under clear-sky conditions (Foga et al., 2017). Therefore, it is important to select proper confidence scores to retain as many clear-sky sea ice pixels as possible while minimizing the number of FN cases. In this study, pixels with medium and high confidence scores for clouds and for cirrus, respectively, were discarded prior to Landsat-8 SIC production to avoid cloud and cirrus contamination. In addition, as suggested in Foga et al. (2017), dilated cloud pixels, which are clear pixels completely surrounded by cloud pixels, were also masked to prevent contamination by cloud edges where cloud detection uncertainty is high.

It is important to note that CFMask over ice-covered surface types has lower accuracy than other surface types (Foga et al., 2017; Qiu et al., 2019). Therefore, an additional step for cloud-masking quality assessment is designated in this study. In this step, a visual inspection was performed by comparing the cloud mask array, which is constructed by masking cloud, cirrus, cloud shadow, and dilated cloud pixels, from each Landsat-8 image with the corresponding true-colour image to identify the cases of FN, false-positive (FP; clear pixels that are mistaken as cloud pixels), true-negative (TN; clear pixels correctly detected as clear pixels), and true-positive (TP; cloud pixels correctly detected as cloud pixels) pixels in the Landsat-8 image. From this additional step, Landsat-8 images were sorted into four categories depending on the assessed quality of cloud masking. Images with the existence of FN cloud pixels in the cloud mask array, which indicate underestimated cloud cover, were labelled as category 1 (C1). Images dominated by FP cloud pixels, which occur in cases of overestimated cloud cover, were tagged as C2. Images dominated by TP cloud pixels, which correspond to correctly estimated cloud cover for cloudy-sky conditions, were labelled as C3. Images dominated by TN cloud pixels, which correctly estimate clear sky, were labelled as C4. For images under C2 (i.e. overestimated cloud coverages with medium confidence scores for clouds and high confidence scores for cirrus), the cloud mask array was regenerated with a higher confidence score (high-confidence clouds and cirrus) and was visually inspected against the true-colour image to determine the adequacy of the higher-confidence-score cloud mask as follows: if any FN cloud pixels were present in the higher-confidence cloud mask, the original confidence score (i.e. medium for clouds and high for cirrus) was used to mask the clouds. Further details of the visual screening step are provided in Appendix B.

In this study, for Landsat-8 images that were labelled as C2, C3, and C4, Landsat-8 pixels that remain after the application of CFMask were assumed to be clear-sky pixels (i.e. “clear-pixel assumption”). However, for Landsat-8 images labelled as C1, the clear-pixel assumption is not valid because the C1 category underestimates clouds through CFMask according to the visual inspection step, which implies that the associated error due to the underestimated cloud cover in SIC calculation is expected. Therefore, possible error from the presence of unmasked cloud pixels in C1 is further evaluated in Sect. 5.1. The number of Landsat-8 images under the four categories over the 12 regions is provided in Table 2, and the assessed cloud mask quality (i.e. C1, C2, C3, and C4) for each Landsat-8 image is provided in the variable under the name “cloud_contamination_category” in the produced Landsat-8 SIC dataset in order to allow for quality-controlling of the data in its usage.

Table 2 The number of Landsat-8 images for the four cloud mask categories (i.e. C1: underestimated cloud cover; C2: overestimated cloud cover; C3: correctly estimated cloud cover for cloudy sky; and C4: correctly estimated cloud cover for clear sky) over the 12 regions of the Arctic Ocean during the periods of January 2020–December 2022.

^a Underestimated cloud cover. ^b Overestimated cloud cover. ^c Correctly estimated cloud cover for cloudy sky. ^d Correctly estimated cloud cover for clear sky.

3.2 Ice and open-water classification

Classification of a Landsat-8 pixel as ice or open water was performed by applying thresholds to the top-of-atmosphere (TOA) reflectivity at band 5 (NIR) and the normalized difference snow index (NDSI). To do this, first, the reflectivity of a Landsat-8 pixel, which is stored as a 16 bit digital number in the Landsat-8 Collection 2 Level 1 dataset, was scaled to TOA reflectivity using the following equation (Zanter, 2019):

$$\begin{array}{}\text{(2)}& {\mathit{\rho }}_{\mathit{\lambda }}={\displaystyle \frac{M_{\mathit{\rho }}{Q}_{\mathrm{DN}}+A_{\mathit{\rho }}}{\sin \left({\mathit{\theta }}_{\mathrm{SE}}\right)}},\end{array}$$

where M_ρ and A_ρ are the multiplicative and additive scale factors, θ_SE is the solar elevation angle, and Q_DN is the TOA reflectivity of the Landsat-8 pixel in a 16 bit digital number format.

Then, the NDSI was calculated from the scaled reflectivities as follows:

$$\begin{array}{}\text{(3)}& \mathrm{NDSI}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{5}}-{\mathit{\rho }}_{\mathrm{6}}}{{\mathit{\rho }}_{\mathrm{5}}+{\mathit{\rho }}_{\mathrm{6}}}},\end{array}$$

where ρ₅ and ρ₆ are the TOA reflectivities at bands 5 (NIR) and 6 (SWIR) of the OLI sensor, respectively.

The steps for differentiating ice and open-water pixels and for removing possible cloud pixels are shown in Fig. 4. The first entails using the ρ₅ criterion to detect open-water pixels, which have lower reflectivity at band 5 compared to that over ice or cloud pixels. The next step is using the NDSI criterion in order to detect ice pixels, which have a higher NDSI than cloud pixels due to the higher reflectivity of ice at band 5 and the lower reflectivity of ice at band 6 compared to the cloud reflectivities (Hall et al., 1995; Riggs et al., 1996, 1999). The NDSI criterion for the separation of ice and cloud pixels was kept in order to reinforce the cloud removal process in addition to CFMask, as explained in Sect. 3.1. In this study, the thresholds for ρ₅ and NDSI were selected to be 0.08 and 0.45, respectively (Liu et al., 2016; Tanaka and Lu, 2023).

As mentioned in Sect. 3.1, the clear-pixel assumption was applied during the classification of Landsat-8 images labelled as C2, C3, and C4. Accordingly, the performance of the classification based on ρ₅ and NDSI with the selected thresholds was evaluated over clear-sky pixels using the surface classification data from Landsat-8 images (Kern, 2021) mentioned in Sect. 2.3 as reference data. The values of ρ₅ and NDSI were collected separately over open-water and ice pixels in the reference data, and classification over the collected pixels was performed following the procedure in Fig. 4. From the distributions of ice and open-water pixels in the two-dimensional histogram between NDSI (y axis) and ρ₅ (x axis) in Fig. C1, it can be seen that ice and open water are well differentiated by the selected threshold values of ρ₅ and NDSI, respectively (Fig. C1). In addition, for quantitative assessment of the performance of ice and open-water classification, the recall was computed for the open-water and ice categories using the classification result summarized in Table 3 and Eq. (4).

$$\begin{array}{}\text{(4)}& {\mathrm{RC}}_{X\text{as}X}={\displaystyle \frac{N_{X\text{as}X}}{N_{X\text{as}X}+N_{X\text{as}\sim X}}}\end{array}$$

In the above, N_X as X and N_{X as ∼X} are the number of pixels in category X classified as X (TP) and the number of pixels in category X classified as not X (FN), respectively. With the designated thresholds, the recall was found to be 98.94 % for water and 97.67 % for ice. FN classifications of ice into open water can cause negatively biased SIC. The bias due to such classification error was estimated to be 2.33 % from the percentage of the number of ice pixels that were classified as open water in Table 3. Conversely, FN classification of open water into ice can cause positively biased SIC, which was estimated to be 1.06 % from the value in Table 3. Misclassification of ice or open-water pixels into cloud pixels from the application of the NDSI threshold rarely occurred for both the ice and open-water categories. Thus, it can be concluded that the classification method used in this study is accurate over clear-sky pixels. Furthermore, the error from ice and water classification under clear-sky conditions is within the uncertainty range of Landsat-8 SIC, which is discussed in Sect. 4.3.

This classification result may not be applicable for Landsat-8 images tagged as C1 (i.e. underestimated cloud cover) because, as mentioned, in Sect. 3.1, such images do not consist solely of clear-sky pixels but contain cloud pixels undetected by CFMask. Therefore, for Landsat-8 images labelled as C1, in order to understand possible errors in SIC calculation from the designated classification method, it is necessary to evaluate the performance of classification over the undetected cloud pixels. This is discussed further in Sect. 5.1.

Figure 4 The process for separating ice, open water, and possible unmasked clouds using the ρ₅ and NDSI criteria, where ρ₅ and NDSI are the TOA reflectivity at band 5 of the OLI sensor and the normalized difference snow index, respectively.

Table 3 The number of classified pixels for open water, ice, and cloud from the suggested method and surface classification reference data (Kern, 2021). The original categories in the reference data are shown in the rows, and the classified categories from the method in Fig. 4 are shown in the columns. The values inside the parentheses indicate the percentage of pixels from the original category that are classified into open water, ice, and cloud.

3.3 Projection and computation of SIC

After the ice and open-water classification for the selected Landsat-8 pixels, the classified pixels were projected onto the target grid system of the NSIDC Polar Stereographic grid with 6.25 km resolution. The number of ice and open-water pixels within each 6.25 km × 6.25 km grid cell was used to compute SIC for the grid cell according to

$$\begin{array}{}\text{(5)}& \mathrm{SIC}={\displaystyle \frac{N_{\mathrm{ice}}}{N_{\mathrm{ice}}+N_{\mathrm{water}}}}\times \mathrm{100}\left(\mathit{\%}\right),\end{array}$$

where N_ice and N_water are the number of Landsat-8 pixels classified as, respectively, ice and water within each 6.25 km × 6.25 km grid cell. It is noted that some of the grid cells with 6.25 km resolution are not entirely filled by Landsat-8 pixels at the edges of a Landsat-8 image and/or near cloud-masked regions. In this study, this kind of grid cell is referred to as a “partially covered grid cell”. Therefore, SIC computed in such a grid cell is unlikely to be representative of the actual ice coverage over the area covered by the grid cell. To avoid this caveat, a minimum threshold in the number of Landsat-8 pixels for a single 6.25 km × 6.25 km grid cell (N_critical) was applied prior to the computation of SIC. In this study, a specific value of N_critical was introduced as the minimum threshold, which is discussed in the following subsection.

3.4 Sensitivity test and uncertainty analysis

The sensitivity of Landsat-8 SIC to the prescribed thresholds of ρ₅ and NDSI was investigated for each cloud contamination category. In doing so, for each of the four cloud contamination categories (i.e. C1, C2, C3, and C4), 10 scenes were randomly sampled over all 12 regions (Fig. 2); thus, 120 scenes were used for each cloud contamination category for the sensitivity test (see Table S5 in the Supplement for the full list of scenes used for the sensitivity test). SIC values over the selected scenes were calculated using Eq. (5) based on classification results with NDSI and ρ₅ thresholds perturbed by their uncertainties. Values of 0.015 and 0.016 were assigned as the uncertainties of ρ₅ and ρ₆, respectively, following Pinto et al. (2020), which provides the root mean squared differences of the Landsat-8 TOA reflectivities and in situ observed reflectivities at bands 5 and 6. The uncertainty of NDSI was calculated using the Gaussian error propagation method, which can be written for NDSI as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{\mathrm{NDSI}}^{\mathrm{2}}={\left({\displaystyle \frac{\partial \mathrm{NDSI}}{\partial {\mathit{\rho }}_{\mathrm{5}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{5}}}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathrm{NDSI}}{\partial {\mathit{\rho }}_{\mathrm{6}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{6}}}^{\mathrm{2}},\end{array}$$

where ${\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{5}}}$ and ${\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{6}}}$ are the uncertainties of Landsat-8 TOA reflectivities at bands 5 and 6, respectively. Substituting Eq. (3) for NDSI in Eq. (6), the analytical form of the uncertainty in NDSI can be expressed as the following:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{NDSI}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{4}{\mathit{\rho }}_{\mathrm{6}}^{\mathrm{2}}}{({\mathit{\rho }}_{\mathrm{5}}+{\mathit{\rho }}_{\mathrm{6}}{)}^{\mathrm{4}}}}{\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{5}}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{4}{\mathit{\rho }}_{\mathrm{5}}^{\mathrm{2}}}{({\mathit{\rho }}_{\mathrm{5}}+{\mathit{\rho }}_{\mathrm{6}}{)}^{\mathrm{4}}}}{\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{6}}}^{\mathrm{2}}.\end{array}$$

From Eq. (7), with ${\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{5}}}=\mathrm{0.015}$ and ${\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{6}}}=\mathrm{0.016}$, a value of 0.05 was assigned as the uncertainty of NDSI, which is the median value of σ_NDSI computed over 480 randomly selected Landsat-8 scenes. For the four cloud contamination categories, mean values of SIC calculated with the perturbed thresholds of 0.45 ± 0.05 and 0.08 ± 0.015 for the NDSI and ρ₅, respectively, are provided in Fig. 5. With the perturbation of  ± 0.015 for the ρ₅ threshold, mean SIC values from C1, C2, C3, and C4 vary by ∓0.641 %, ∓0.495 %, ∓0.665 %, and ∓0.402 %, respectively (blue lines in Fig. 5). With a perturbation of ±0.05 for the NDSI threshold, mean SIC values from C1, C2, C3, and C4 varied by ∓0.111 %, ∓0.002 %, ∓0.007 %, and ∓0.002 %, respectively (red lines in Fig. 5). The calculated SIC values are more sensitive to the ρ₅ threshold relative to the NDSI threshold because the ρ₅ threshold is responsible for separating open water and ice. It is noted that the sensitivity of SIC values to the NDSI threshold is 2 orders of magnitude higher for scenes labelled as C1 than for C2, C3, and C4. The very low sensitivity of SIC values to the NDSI threshold for scenes labelled as C2, C3, and C4 leads us to infer that cloud pixels in such scenes had been successfully masked by CFMask prior to the ice or water classification described in Sect. 3.2. However, the relatively higher sensitivity of SIC values to the NDSI threshold for scenes labelled as C1 leads us to infer that undetected cloud pixels remained after the application of CFMask and that such cloud pixels had been further removed by the NDSI threshold.

Gaussian error propagation was also used to estimate the uncertainty of Landsat-8 SIC according to the following:

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{\mathrm{SIC}}^{\mathrm{2}}={\left({\displaystyle \frac{\partial \mathrm{SIC}}{\partial \mathrm{NDSI}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{NDSI}}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathrm{SIC}}{\partial {\mathit{\rho }}_{\mathrm{5}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\rho }}_{\mathrm{5}}}^{\mathrm{2}}\left(\mathit{\%}\right),\end{array}$$

where σ_x and $\frac{\partial \mathrm{SIC}}{\partial x}$ are the uncertainty of x and the sensitivity of SIC to x, respectively. The sensitivities for the two variables (i.e. ρ₅ and NDSI) were computed numerically from the mean SIC variation observed in the sensitivity test (see Tables S6, S7, S8, and S9 in the Supplement for the computed values of sensitivity). In addition, in order to check the relative contribution of each variable to the overall uncertainty in SIC, a contribution factor (CF_x) was defined and calculated for the two variables as follows:

$$\begin{array}{}\text{(9)}& {\mathrm{CF}}_x={\displaystyle \frac{{\left(\frac{\partial \mathrm{SIC}}{\partial x}\right)}^{\mathrm{2}}{\mathit{\sigma }}_x^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{SIC}}^{\mathrm{2}}}}\times \mathrm{100}\left(\mathit{\%}\right).\end{array}$$

The estimated uncertainty of Landsat-8 SIC (σ_SIC) produced in this study was less than 1 %, on average, for all four cloud contamination categories, and the ρ₅ threshold contributes to about 99 % of the uncertainty for C2, C3, and C4 and to about 97 % of the uncertainty for C1 in the SIC calculation. Further discussion of the uncertainty of Landsat-8 SIC is presented in Sect. 4.3.

As mentioned in Sect. 3.3, SIC computed from partially covered grid cells may not be representative of actual ice coverage over the entire grid cell, and the corresponding uncertainty of SIC estimates in such grid cell can be as large as the fraction of the uncovered areas. In order to circumvent such a problem, in this study, N_critical was determined to be 0.99 × N_max, where N_max is the maximum number of Landsat-8 pixels within a 6.25 km × 6.25 km grid cell.

Figure 5 Mean of Landsat-8 SIC values for (a) C1 (i.e. underestimated cloud cover), (b) C2 (i.e. overestimated cloud cover), (c) C3 (i.e. correctly estimated cloud cover for cloudy sky), and (d) C4 (i.e. correctly estimated cloud cover for clear sky) derived from the selected scenes under perturbed thresholds for NDSI (red) and ρ₅ (blue), where ρ₅ and NDSI are the TOA reflectivity at band 5 of the OLI sensor and the normalized difference snow index, respectively.

3.5 Application of land and region masks

In order to circumvent potential contamination of land signals, in this study, SIC pixels generated over non-ocean regions were masked using the surface mask described in Sect. 2.2. The region mask was applied in addition to the surface mask to obtain SIC products catalogued into the 12 regions. If all SIC pixels in a Landsat-8 scene were masked by the combination of land, region, and cloud masks, the scene was removed from the SIC dataset.

4 Result

4.1 Landsat-8 SIC dataset

Out of 15 286 Landsat-8 level-1 images collected over pan-Arctic areas during the study period, the numbers of Landsat-8 images used for the calculation of SIC values for the categories of C1, C2, C3, and C4 were 6336 (41.4 %), 549 (3.6 %), 4389 (28.1 %), and 3123 (20.4 %), respectively. The remaining 989 images (6.5 %) were removed from the combination of surface, region, and cloud masks. For each of the 12 regions, the number of Landsat-8 scenes generated into Landsat-8 SIC (N_scene) and the number of produced Landsat-8 SIC pixels (N_pixel) for each cloud contamination category during the study period are shown in Fig. 6, along with the mean and standard deviation of SIC (see Table S10 in the Supplement for values). The total number of Landsat-8 SIC pixels produced during the study period was 2 934 399.

Figure 6 (a) The mean SIC, (b) the standard deviation of SIC, (c) the number of Landsat-8 scenes used for SIC production (N_scene), and (d) the number of Landsat-8 SIC pixels produced (N_pixel) over the 12 regions. The black, red, blue, and green bars indicate values for categories C1 (i.e. underestimated cloud cover), C2 (i.e. overestimated cloud cover), C3 (i.e. correctly estimated cloud cover for cloudy sky), and C4 (i.e. correctly estimated cloud cover for clear sky), respectively.

4.2 Qualitative evaluation for Landsat-8 SIC under four cloud contamination categories

Figure 7 shows the Landsat-8 true-colour image (first column of Fig. 7); classification map of ice, open water, and removed cloud pixels (second column of Fig. 7); and Landsat-8 SIC at 6.25 km resolution (third column of Fig. 7) for the four selected cases. Ice and open-water pixels, which were differentiated following the methods explained in Sect. 3.2, are shown as the white and blue pixels, respectively. Cloud pixels removed in both CFMask and the NDSI criterion are shown as the cyan pixels. Cloud pixels removed from CFMask but undetected using the NDSI criterion are shown as the purple pixels. Cloud pixels removed from the NDSI criterion but undetected in CFMask are shown as the red pixels. SIC values were only estimated over grid cells that satisfy N>N_critical; therefore, grid cells with more than 1 % of their area covered with cloud pixels or grid cells located near the edges of a Landsat-8 scene have no SIC values. In addition, areas close to the coastline (within 6.25 km) are masked in the Landsat-8 SIC maps presented in Fig. 7.

The first case is an example of the underestimated cloud cover (i.e. C1) on 13 March 2022 over the Kara Sea (first row of Fig. 7), where cloud pixels observed in the lower-left area of Fig. 7a were not removed by CFMask (cyan and purple pixels in Fig. 7b). However, for this particular scene, most of such undetected cloud pixels were removed from the application of the NDSI criterion (red pixels in Fig. 7b), and, therefore, the produced SIC was estimated only over clear-sky area (Fig. 7c). The second case is an example of the overestimated cloud cover (i.e. C2) on 17 March 2021 over the Barents Sea (second row of Fig. 7), where FP cloud pixels are densely distributed in the upper-left area of Fig. 7e. It is shown that SIC values were not estimated for grid cells with such wrongly masked pixels (Fig. 7f). The third is an example of correctly estimated cloud cover for cloudy sky (i.e. C3) on 26 June 2022 over the Kara Sea (third row of Fig. 7), where the position of cloud pixels removed from CFMask (cyan and purple pixels in Fig. 7h) coincides well with the location of cloud presented in the true-colour image (Fig. 7g). The fourth case is an example of correctly estimated cloud cover for clear sky (i.e. C4) on 15 June 2022 over the Beaufort Sea (fourth row of Fig. 7), where no clouds are observed in both the true-colour image (Fig. 7j) and the classification map (Fig. 7k).

For all four cases, over clear-sky pixels, discrimination between open-water pixels (blue pixels in Fig. 7b, e, h, and k) and ice pixels (white pixels in Fig. 7b, e, h, and k) based on the ρ₅ thresholds coincided well with the locations of open water and ice observed from the true-colour images (first column in Fig. 7). Therefore, it can be concluded that the ice–water classification in this study is successfully done and that the calculated SIC values correspond well to the classification results (third column in Fig. 7). In addition, cloud pixels only detected from the NDSI criterion (red pixels in second column in Fig. 7) are rarely observable for the cases of C2, C3, and C4, which further demonstrates the validity of the assumption that all cloud pixels had been removed prior to ice–water classification in Sect. 3.2 for the Landsat-8 scenes under the three categories.

Figure 7 Example of (a, d, g, j) original Landsat-8 true-colour image; (b, e, h, k) classification map of ice (white), open water (blue), and cloud (cyan, purple, and red); and (c, f, i, l) Landsat-8 SIC values at 6.25 km resolution on (first row) 22 March 2022 over the Kara Sea, (second row) 17 March 2021 over the Barents Sea, (third row) 26 June 2022 over the Kara Sea, and (fourth row) 15 June 2022 over the Beaufort Sea. From top to bottom, the select cases correspond to the cloud contamination categories 1, 2, 3, and 4 respectively. SIC values are not estimated over areas of cloud mask (cyan, purple, and red pixels in the middle column), and SIC values near the coastal area (6.25 km) are masked in this figure. The true-colour images in (a, d, g, j) are courtesy of the US Geological Survey. © Authors for panels (b), (c), (e), (f), (h), (i), (k), and (l). Distributed under the Creative Commons Attribution 4.0 License. Landsat-8 for panels (a), (d), (g), and (j).

4.3 Uncertainty of Landsat-8 SIC

The estimated σ_SIC from all 480 selected scenes in Sect. 3.4 was less than 1 %, and the ρ₅ threshold was found to be responsible for more than 99 % of σ_SIC. The uncertainties (i.e. σ_SIC) estimated separately for different regions or different cloud contamination categories all remained within 1 % and had similar contribution ratios, with the ρ₅ threshold being the dominant factor contributing to σ_SIC (see Tables D1, D2, and D3 in Appendix D for the exact values). Thus, σ_SIC seems to be independent of region or cloud contamination label. However, σ_SIC was found to be dependent on the SIC value itself. Figure 8 shows the variation in σ_SIC with respect to the SIC range, illustrating that the lowest uncertainty is ∼ 0.2 % in the SIC values ranging from 0 % to 10 % and from 90 % to 100 %, while the highest uncertainty of 4.5 % is observed in SIC values ranging from 50 % to 60 % (see Table D4 for exact values). Still, the ρ₅ threshold explains most of the uncertainty regardless of the SIC values. In spite of the relatively high uncertainty in Landsat-8 SIC between 20 % and 80 %, the product can still be used for validation purposes because most PMW SIC products exhibit much larger uncertainties over such SIC ranges of up to ±12 % in the winter (Ivanova et al., 2014) and ±20 % in the summer (Meier and Notz, 2010).

Figure 8 (a) Uncertainties in Landsat-8 SIC values (σ_SIC) and (b) contributions of the ρ₅ (blue) and the NDSI thresholds (red) to the estimated uncertainties for different SIC ranges, where ρ₅ and NDSI are the TOA reflectivity at band 5 of the OLI sensor and the normalized difference snow index, respectively. Dark- and light-coloured bars indicate the uncertainty and contributions computed from all 480 scenes and separately for each SIC sub-range, respectively.

5 Discussion

5.1 Possible errors in SIC produced from Landsat-8 images labelled as C1

As mentioned in Sect. 3.1 and 3.2, the clear-pixel assumption, which assumes that all cloud pixels in a Landsat-8 image have been removed by the application of CFMask, is not valid for Landsat-8 images labelled as C1 in Sect. 3.1. Therefore, for Landsat-8 SIC associated with the C1 category, it is necessary to evaluate the possible uncertainty in SIC induced by unremoved cloud pixels (i.e. underestimated cloud cover). The evaluation was performed as follows: from Landsat-8 images under the C1 category, those having 100 % cloud cover based on visual inspection but less than 100 % cloud cover based on CFMask were selected. From these images, the ρ₅ and NDSI values were collected from pixels that were not masked by CFMask (i.e. undetected cloud pixels). Classification following the process illustrated in Fig. 4 was performed over the collected undetected cloud pixels to quantitatively assess the possible errors in SIC estimated over such pixels. A total of 6 721 605 undetected cloud pixels were used in this evaluation, and the name and location of the Landsat-8 images used are shown in Fig. S3 and in Table S11 in the Supplement.

The classification result is summarized in Table 4 and Fig. C2. From the distribution of the unmasked cloud pixels in the two-dimensional histogram between NDSI (y axis) and ρ₅ (x axis) in Fig. C2, it can be seen that the NDSI criterion used in this study reinforces the cloud removal process by filtering cloud pixels that were undetected by CFMask. However, even with the additional procedure to remove remaining cloud signals (i.e. the NDSI criterion), 8.54 % of the undetected cloud pixels are still classified as open water and/or ice. It is noted that such cloud pixels (i.e. cloud pixels undetected by both CFMask and the NDSI criterion) were predominantly classified as ice (Table 4). Therefore, it can be inferred that the undetected cloud pixels in a Landsat-8 image can induce positively biased SIC, and, thus, for SIC values produced from Landsat-8 images labelled as C1 over which the clear-pixel assumption is invalid, the error from ice–water classification is estimated to be large as 8.54 % from the percentage of cloud pixels classified as ice in Table 4. The possibility of such large uncertainties should be taken into account when using Landsat-8 SIC labelled as C1.

Table 4 The number of cloud pixels that were undetected using the C Function of Mask (CFMask) classified into open water, ice, and cloud from the application of the procedure in Fig. 4. The scenes used for the evaluation belong to the C1 (i.e. underestimated cloud cover) category from the method described in Sect. 3.1. The values inside the parentheses indicates the percentage of pixels that belong to each category.

5.2 Evaluation of Landsat-8 SIC using ice charts

The accuracy of Landsat-8 SIC produced in this study was evaluated using ice charts provided by MET Norway as reference. For quantitative comparison, ice charts were collocated into the grid system of Landsat-8 SIC (i.e. PSR grid with 6.25 km resolution) as follows: data points on the ice chart within each 6.25 km × 6.25 km grid cell were collected, and the SIC mean value from the collected data points were taken as the representative SIC value of the ice chart for the 6.25 km × 6.25 km grid cell. It is important to note that SIC values in the original ice charts are not normally defined SIC values in satellite remote sensing but contain uncertainties represented by the ice concentration range defined in Table A1. Therefore, it is necessary to consider the propagation of uncertainty in the collocation process. Uncertainty of the collocated ice chart SIC was estimated by taking the average of the uncertainty in ice chart data points collected from each 6.25 km × 6.25 km grid cell. To avoid the influence of land contamination, a 6.25 km coastal area mask was applied to both SIC values prior to the comparison. The number of collocated data points used in the evaluation was 45 547.

From Fig. 9, a good linear relationship (i.e. correlation coefficient of 0.96) between Landsat-8 SIC and ice chart SIC is observed. The spread (i.e. 20th and 80th percentiles) of Landsat-8 SIC for ice chart SIC sub-ranges, which are shown as vertical red lines in Fig. 9a, was larger in SIC ranging from 20 % to 80 % relative to other ranges, which is likely a consequence of the wider concentration intervals assigned to the 20 %–80 % SIC values in the original ice chart (Table A1). In addition, SIC from the ice charts tends to be higher than that found from Landsat-8, and the bias is more pronounced in the lower SIC range. This type of state-dependent overestimation of SIC from ice charts has been reported in the previous works of Tonboe et al. (2016) and Cheng et al. (2020), showing that overestimation of SIC from ice charts is largest in the lower SIC range due to the “better-safe-than-sorry” practices of the ice-charting community. For quantitative comparison of the bin-wise mean biases in Landsat-8 SIC relative to ice chart SIC, bin-averaged SIC values from Landsat-8 (red triangle in Fig. 9b) and from ice charts (blue circles in Fig. 9b) were plotted, along with their respective uncertainties. Uncertainties in Landsat-8 SIC values over the SIC sub-ranges were taken as the values from Table D4. Except for the 70 %–80 % SIC interval, Landsat-8 SIC values were negatively biased in relation to ice chart SIC values. However, the mean biases for all SIC sub-ranges were found to be within the uncertainty ranges estimated for each product.

Figure 9 (a) Scatterplot of bin-wise mean Landsat-8 SIC values and ice chart SIC sub-ranges. The bin-wise mean SIC values are shown as red triangles, and the 20th and 80th percentiles are shown as the vertical red lines. The values for the number of data points (N), root-mean-square error (RMSE), bias, and Pearson correlation coefficient (R) are presented. (b) For the same SIC intervals as (a), the bin-wise mean SIC values of Landsat-8 (red triangle) and ice charts (blue circle) and their respective uncertainties (vertical lines). The uncertainties of Landsat-8 SIC are taken from the values in Table D4.

5.3 Evaluation of Landsat-8 SIC over melt ponds

The evaluation of variation in Landsat-8 SIC values due to melt pond presence was performed using the MPF dataset (Niehaus and Spreen, 2022) described in Sect. 2.5 as reference data for melt ponds. As mentioned in Sect. 2.5, a total of six Landsat-8 scenes obtained from the periods of July 2020, August 2020, and July 2021 were used in the evaluation. The evaluation was conducted as follows: first, the collocation of the MPF dataset into the grid system of Landsat-8 SIC was performed. This was done by collecting MPF data points within each 6.25 km × 6.25 km grid cell and taking the mean value of the collected MPF data points as the MPF values for each corresponding grid cell at 6.25 km resolution. In addition, from the OW masks in the MPF dataset, SIC values (SIC_MPF) were computed in the grid system of Landsat-8 SIC following the same method introduced in Sect. 3.3. Second, in order to remove the effects of SIC variation from the evaluation, the corresponding Landsat-8 SIC and MPF data points were collected when the data points satisfied SIC_MPF=100 %. The number of collected data points is 98. From the collected data points, the net ice surface fraction (C_net) was computed as follows:

$$\begin{array}{}\text{(10)}& C_{\mathrm{net}}=(\mathrm{1}-\mathrm{MPF})\times {\mathrm{SIC}}_{\mathrm{MPF}},\end{array}$$

where MPF is the melt pond fraction, and SIC_MPF is the estimated SIC value from the MPF dataset. Since SIC_MPF was fixed at 100 %, in this study, the variation in C_net can be considered to be driven solely by the variation in MPF.

The robustness of the Landsat-8 SIC in relation to the presence of melt ponds is illustrated in Fig. 10, which is a scatterplot of the collected C_net (x axis) and Landsat-8 SIC (y axis). In this plot, the MPF varied from 0 % to 33 %, and, therefore, the computed values of C_net are ranged from 67 % to 100 %. However, SIC values estimated from Landsat-8 are observed to be nearly independent of the varying C_net (statistically insignificant correlation coefficient of 0.11) and, thus, nearly independent of MPF. Although a few Landsat-8 SIC values are observed to be affected by melt pond presence (data points highlighted in red in Fig. 10), which can be expected because melt ponds are not easily distinguished from open water, the number of such data points is very small (only 4 data points out of 98 data points). It is noted that, on average, the deviation from 100 % ice concentration computed from the data points shown in Fig. 10 was less than 1 %. Therefore, it can be inferred that the impact of melt pond presence in SIC calculation using Landsat-8 imagery is small and that the proposed algorithm for SIC production in this study is robust regardless of surface melting.

Figure 10 Scatterplot of net ice surface fraction (x axis) and Landsat-8 SIC (y axis). The data points shown satisfy SIC_MPF=100 % and have MPF values that vary from 0 % to 33 %. Data points with Landsat-8 SIC values that deviate by more than 4 % from 100 % ice concentration are highlighted in red. The values for the number of data points (N), the Pearson correlation coefficient (R), and the p value for the correlation coefficient are presented.

5.4 Possible applications of Landsat-8 SIC for assessing PMW-based SIC values

Landsat-8 SIC produced from this study can be utilized to assess the PMW-based SIC values. This section provides the possible application of the constructed Landsat-8 SIC for examining PMW-based SIC values. To do this, Landsat-8 SIC was downscaled to 25 km resolution and compared against SIC values estimated from BT and NT algorithms, both provided in the PSR grid with 25 km resolution and obtained from NSIDC (Meier et al., 2021), for the two selected cases of Landsat-8 scenes acquired during melting (21 July 2022 over the Laptev Sea) and freezing (4 March 2022 over the Chukchi Sea) seasons, respectively. To avoid the influence of land contamination, a coastal area mask, which was also downscaled to 25 km resolution, was applied before the comparison.

Figure 11 illustrates the spatial distributions of the three different SIC values, differences in SIC values of BT and NT from Landsat-8, and scatterplots of BT and NT SIC against Landsat-8. For the case of the melting season (top two rows in Fig. 11), BT SIC showed a positive bias of 8.95 %, an RMSE of 16.30 %, and a correlation coefficient of 0.92 (Fig. 11d) in relation to Landsat-8 SIC, while SIC retrieved from the NT algorithm is negatively biased in relation to Landsat-8 SIC by −5.21 %, with an RMSE of 14.35 % and a correlation coefficient of 0.94 (Fig. 11h). It is interesting to note that BT SIC values are positively (negatively) biased in relation to Landsat-8 SIC for lower (higher) concentrated ice areas (Fig. 11c), while opposite patterns are observed for NT SIC values (Fig. 11g). Both PMW-based SIC values show the largest disagreement with Landsat-8 SIC near the edges of pack ice (i.e. boundaries between sea ice and open water).

For the scene in the freezing season (bottom two rows in Fig. 11), the BT and NT algorithms produced nearly 100 % SIC values for all grids in this case, while Landsat-8 shows lower SIC values in regions coinciding with the leads in the pack ice observed from the true-colour image. As a result, positive biases were observed near the position of the lead (Fig. 11k and o), and mean biases for the BT and NT algorithms were 0.83 % and 0.53 %, respectively. RMSEs of BT and NT SIC were calculated to be 1.35 % and 0.81 %, respectively, which are lower than the RMSEs evaluated during the melting season for the two SIC algorithms (Fig. 11i and p).

Figure 11 Geographical distributions of (a, i) original Landsat-8 true-colour image, (e, m) Landsat-8 SIC, (b, j) SIC from the bootstrap (BT) algorithm, (f, n) SIC from the NASA Team (NT) algorithm, (c, k) difference in SIC values between BT and Landsat-8, (g, o) difference in SIC values between NT and Landsat-8, and scatterplot (d, l) between Landsat-8 SIC and SIC from BT and (h, p) between Landsat-8 SIC and SIC from NT. The values of root-mean-square error (RMSE), bias, and Pearson correlation coefficient (R) are presented with the scatterplots. Upper two panels for 21 July 2022 (melting season) over the Laptev Sea and for 4 March 2022 over the Chukchi Sea. The true-colour images are courtesy of the US Geological Survey. The SIC retrievals from the BT and NT algorithms were obtained from Meier et al. (2021).

Figure 12 Variables in the Landsat-8 SIC NetCDF. The scene is from 12 June 2021 over the Canadian Archipelago. For the (d) sub-region mask, “in” and “out” denote grid cells located inside and outside the designated region, respectively.

Table 5 Variables in the Landsat-8 SIC NetCDF file, the name of the variables, and the flag values for each variable.

6 Code and data availability

The Landsat-8 SIC dataset can be downloaded at https://doi.org/10.5281/zenodo.10973297 (Jung et al., 2024). Datasets generated for each Arctic region can be found in sic_landsat08_{region name}.nc. The datasets are stored in NetCDF format and can be accessed using software including Python, MATLAB, and QGIS. Along with the SIC values, the N, coastal mask, and region mask are also provided in a 1792 × 1216 array format. The cloud contamination categories and names of the original Landsat-8 files are also provided for each scene. Variables in the NetCDF file are visualized in Fig. 12. Fill values were assigned to grids outside the coverage of a Landsat-8 scene, grids over land, or grids masked by clouds (black grids in Fig. 12a, b). Descriptions of each variable and the fill and flag values are summarized in Table 5.

Datasets used to produce and validate the Landsat-8 SIC are listed as follows:

• The Landsat-8 Collection 2 Level 1 product and the corresponding true-colour images are accessible from the United States Geological Survey Earth Explorer at https://earthexplorer.usgs.gov/ (last access: 16 March 2024; DOI: https://doi.org/10.5066/P975CC9B, Earth Resources Observation and Science (EROS) Center, 2020).

• The Arctic and Antarctic Regional Masks for Sea Ice and Related Data Products version 1, used to mask non-ocean areas and to distinguish regions, can be accessed at https://doi.org/10.5067/CYW3O8ZUNIWC (Meier and Stewart, 2023a).

• Landsat surface type over water from supervised classification of surface broadband albedo estimates (version_2021_fv0.01), used to test the performance of the ice and open-water classification, can be accessed at https://doi.org/10.25592/uhhfdm.9181 (Kern, 2021).

• The Arctic Ocean – Sea Ice Concentration Charts – Svalbard and Greenland ice charts used to evaluate the produced Landsat-8 SIC can be accessed at https://doi.org/10.48670/moi-00128 (E.U. Copernicus Marine Service Information (CMEMS), 2024).

• Melt pond fraction on Arctic sea-ice from Sentinel-2 satellite optical imagery (2017–2021) used to test the robustness of Landsat-8 SIC over melt ponds can be accessed at https://doi.org/10.1594/PANGAEA.950885 (Niehaus and Spreen, 2022).

• NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration version 4, used to illustrate possible applications of the Landsat-8 SIC dataset, can be accessed at https://doi.org/10.7265/efmz-2t65 (Meier et al., 2021).

The Python codes for Landsat-8 SIC production, sensitivity and uncertainty analysis, ice–water classification evaluation, Landsat-8 SIC validation, and figure generation are accessible at https://doi.org/10.5281/zenodo.12754602 (Jung, 2024). Example data to check the functionalities of each Python code are provided with the code repository.

7 Conclusion

In this study, 3 years (2020–2022) of Landsat-8 data were collected and used to generate sea ice concentration (SIC) datasets in the polar stereographic grid at 6.25 km resolution. A total of 14 297 Landsat-8 images were used to calculate 2 934 399 SIC grid points. Each Landsat-8 SIC is catalogued under a NetCDF file named after the 12 regions.

Each Landsat-8 image was labelled as one of four cloud contamination categories (i.e. C1, C2, C3, and C4) according to the overall quality of the cloud mask over the image. The categories are provided in the variable under the name cloud_contamination_category of the Landsat-8 SIC dataset to allow for selection of SIC values calculated without the interference of cloud signals.

The uncertainty of Landsat-8 SIC was estimated to range from 1 % to 4 % based on the Gaussian error propagation method. In addition, to regulate the potential uncertainty that may arise from the use of partially covered grid cells, SIC was only produced for grid cells with over 99 % of their area covered by Landsat-8 pixels. Evaluation of Landsat-8 SIC using SIC from ice charts shows good linear correlation between the two products and also reveals existence of a negative bias in Landsat-8 SIC. However, the bias was found to be within the uncertainty range of the Landsat-8 and ice chart SIC values. In addition, the production method used for Landsat-8 SIC was found to be robust over melt ponds. Thus, the Landsat-8 SIC values produced in this study can be considered to be reliable estimates of SIC.

Comparison of Landsat-8 SIC against SIC retrievals from the NASA Team (NT) and bootstrap (BT) algorithms for two cases reveals an overall negative bias in NT and a positive bias in BT SIC. The spatial distribution of the bias shows that bias in NT and BT SIC may be related to the SIC values, with NT SIC exhibiting a stronger negative bias in the high-SIC regime and with BT SIC showing a stronger positive bias in the low-SIC regime. This suggests that the Landsat-8 SIC can be used as the reference SIC to generate quantitative error statistics of various passive microwave SIC retrievals over different regions, seasons, and SIC values, which can be used to develop an optimal combination of existing SIC algorithms or to provide realistic observation errors to enhance the performance of sea ice data assimilation.

Future works should aim to extend the temporal and spatial coverage of the current Landsat-8 SIC dataset through the addition of Landsat-8 images from the years 2018 and 2019. In addition, given the large number of Landsat-8 SIC data points generated in this study, the obtained SIC values also have the potential to be used to train deep-learning models in order to retrieve optimal SIC estimates over the Arctic.

Appendix A: SIC values in the MET Norway ice chart

Table A1 Concentration class, concentration interval, fixed concentration value, and concentration range of the operational ice chart produced by MET Norway.

Appendix B: Visual inspection for cloud-masking quality control

In this section, a step-by-step description of the process to perform the visual inspection of Landsat-8 scenes is presented. As defined in Sect. 3.1, each pixel in a Landsat-8 scene can be sorted into the following four categories depending on the state of the cloud mask for the pixel: false negative (FN – cloud pixel mistaken as clear pixel), false positive (FP – clear pixel mistaken as cloud pixel), true negative (TN – clear pixel identified as clear pixel), and true positive (TP – cloud pixel identified as cloud pixel). It is noted that the pixels with an FN are used to calculate SIC values, while the pixels with an FP are not, indicating that the presence of FN pixels can directly introduce errors into the calculated SIC value. Therefore, visual inspection was performed very strictly to detect FN pixels.

Figure B1 outlines the steps taken to perform the visual inspection. The descriptions of each step are provided, along with an example case of a Landsat-8 scene that is categorized into C1 during the section-wise inspection stage.

• Step 1. Generate a portable network graphics (PNG) file of the cloud mask (i.e. cloud mask image).

  For each Landsat-8 scene, a false-colour image for each pixel classified as ice (white pixels in Fig. B2b, d), open water (blue pixels in Fig. B2b, d), cloud (grey pixels in Fig. B2b, d), and fill value (black pixels in Fig. B2b, d) is constructed using the OpenCV module in Python. Ice and open-water pixels are differentiated using the method described in Sect. 3.2. Cloud pixels are classified by masking the medium-confidence cloud, high-confidence cirrus, cloud shadow, and dilated cloud pixels identified by the quality assessment band (i.e. the cloud mask array produced by CFMask).

• Step 2. Conduct a comprehensive inspection of the cloud mask quality.

  The cloud mask image generated in step 1 is visually inspected against the true-colour image to identify sections populated with FN, FP, TN, or TP pixels. This is done in the following order: first, if no cloud pixels are observed in both the cloud mask image and the true-colour image (i.e. all pixels in the image are TN pixels), the scene is labelled as C4. Second, if any cluster of FN pixels is observed, the scene is labelled as C1. Third, if any cluster of FP pixels is identified, the scene is labelled as C2 and passes on to step 3. If the clusters of cloud pixels in the cloud mask image correspond with the position of clouds observed in the true-colour image (i.e. TP pixels) well, the scene is labelled as C3 and passes on to Step 4 (Fig. B2a, b).

• Step 3. Conduct a comprehensive inspection of the cloud mask quality of C2.

  For the scenes passed on to this step (i.e. scenes labelled C2 from step 2), the cloud mask image is recreated using a higher-confidence threshold (i.e. high-confidence cloud, high-confidence cirrus, cloud shadow, and dilated cloud pixels) for the quality assessment band. The new cloud mask image is visually inspected against the true-colour image, and if any cluster of FN pixels is observed, the confidence threshold for the quality assessment band is returned to its initial value (i.e. medium-confidence cloud, high-confidence cirrus, cloud shadow, and dilated cloud pixels). If the observed clusters of cloud pixels in the new cloud mask image correspond with the position of clouds observed in the true-colour image well, the higher-confidence threshold is kept, and the scene passes on to step 4.

• Step 4. Conduct a section-wise inspection of cloud mask quality.

  In this step, the identified clusters of TP pixels are inspected in more detail. For each cluster of TP pixels observed, we zoom in (i.e. to about 1000 × 1000 pixels, with the full-size image being approximately 8000 × 8000 pixels) to the section of the cluster to check for the existence of FN pixels. If any FN pixels are found within the cluster, the scene is labelled as C1 (Fig. B2c, d).

An example of how a Landsat-8 scene may be categorized according to the process described in Fig. B1 is presented using the case of a Landsat-8 scene acquired on 25 March 2022 over the Barents Sea (Fig. B2). First, from step 2 (i.e. the comprehensive-inspection step), visual inspection of the cloud mask image (Fig. B2b) against the true-colour image (Fig. B2a) shows that the position of clouds in the cloud mask array generally corresponds with those observed in the true-colour image well. Therefore, in this step, this scene is labelled as C3 and passes on to step 4 as described above. Next, the section-wise inspection of the cloud mask quality is performed by zooming in on the cloud areas. This is illustrated in Fig. B2c and d, which is a zoomed-in image of the area enclosed by the red rectangle in Fig. B2a and b. Inspection of this subsection shows the presence of unmasked cloud shadow pixels, which results in the erroneous classification of ice as open water. Therefore, during this step, the label of this scene is changed to C1.

The visual inspection was done by the author Hee-Sung Jung, and it took approximately 5–10 min to inspect one Landsat-8 scene for cloud cover.

Figure B1 Processing pipeline of the visual inspection step. Each Landsat-8 image is labelled as C1 (i.e. underestimated cloud cover), C2 (i.e. overestimated cloud cover), C3 (i.e. correctly estimated cloud cover for cloudy sky), or C4 (i.e. correctly estimated cloud cover for clear sky) depending on the observed dominance of true-negative (TN – clear pixels identified as clear pixels), false-negative (FN – cloud pixels mistaken as clear pixels), false-positive (FP – clear pixels mistaken as cloud pixels), and true-positive (TP – cloud pixels identified as cloud pixels) pixels.

Figure B2 The case of a Landsat-8 scene classified as C1 (i.e. underestimated cloud cover) from the visual-inspection step. Shown in the panels are (a) fully sized true-colour image, (b) fully sized cloud mask array, (c) true-colour image of the area enclosed by the red rectangle in (a) and (b), and (d) cloud mask array of the area enclosed by the red rectangle in (a) and (b). The blue, white, grey, and black pixels in (b) and (d) are open-water, ice, cloud, and fill value pixels, respectively. The scene is from 25 March 2022, over the Barents Sea. The true-colour image is courtesy of the US Geological Survey.

Appendix C: Validation of ice and open-water classification

Figure C1 Scatterplot of NDSI and ρ₅ for (a) ice and (b) open water. The values for ice and open-water pixels were collected using the ice–water surface classification map (Kern, 2021) as reference data. The thresholds for NDSI and ρ₅ used in this study are shown by the dashed white lines. The colour bars denote the number of pixels.

Figure C2 Scatterplot of NDSI and ρ₅ for cloud pixels that remain unmasked after the application of CFMask. The pixels are acquired from 10 select Landsat-8 images categorized as C1. The thresholds for NDSI and ρ₅ used in this study are shown by the dashed white lines. The colour bars denote the number of pixels.

Appendix D: Uncertainty of Landsat-8 SIC with respect to region, cloud contamination label, and SIC sub-range

Table D1 Contributions of ρ₅ ($C_{{\mathit{\rho }}_{\mathrm{5}}}$) and NDSI (C_NDSI) to the estimated uncertainty of SIC (σ_SIC) and σ_SIC over all 480 scenes.

Table D2 Contributions of ρ₅ ($C_{{\mathit{\rho }}_{\mathrm{5}}}$) and NDSI (C_NDSI) to the estimated uncertainty of SIC (σ_SIC) and σ_SIC for the 12 regions.

Table D3 Contributions of ρ₅ ($C_{{\mathit{\rho }}_{\mathrm{5}}}$) and NDSI (C_NDSI) to the estimated uncertainty of SIC (σ_SIC) and σ_SIC for the four cloud contamination categories.

Table D4 Contributions of ρ₅ ($C_{{\mathit{\rho }}_{\mathrm{5}}}$) and NDSI (C_NDSI) to the estimated uncertainty of SIC (σ_SIC) and σ_SIC for varying SIC sub-range.