Exploratory Mapping of Blue Ice Regions in Antarctica Using Very High-Resolution Satellite Remote Sensing Data

Abstract

Mapping spatiotemporal changes in the distribution of blue ice regions (BIRs) in Antarctica requires repeated, precise, and high-resolution baseline maps of the blue ice extent. This study demonstrated the design and application of a newly-developed semi-automatic method to map BIRs in the Antarctic environment using very high-resolution (VHR) WorldView-2 (WV-2) satellite images. We discussed the potential of VHR satellite data for the mapping of BIRs in the Antarctic environment using a customized normalized-difference blue-ice index (NDBI) method devised using yellow, green, and near-infrared spectral bands of WV-2 data. We compared the viability of the newly developed customized NDBI approach against state-of-the-art target detection (TD), spectral processing (SP) and pixel-wise supervised (PSC) feature extraction (FE) approaches. Four semi-automatic FE approaches (three existing plus one newly developed) consisting of 16 standalone FE methods (12 existing + four customized) were evaluated using an extensive quantitative and comparative assessment for mapping BIRs in the vicinity of Schirmacher Oasis, on the continental Antarctic coastline. The results suggested that the customized NDBI approach gave a superior performance and the highest statistical stability when compared with existing FE techniques. The customized NDBI generally rendered the lowest level of misclassification (average RMSE = 654.48 ± 58.26 m²), followed by TD (average RMSE = 987.81 ± 55.05 m²), SP (average RMSE = 1327.09 ± 127.83 m²) and PSC (average RMSE = 2259.43 ± 115.36 m²) for mapping BIRs. Our results indicated that the use of the customized NDBI approach can greatly improve the semi-automatic mapping of BIRs in the Antarctic environment. This study presents the first refined map of distribution of BIRs around the Schirmacher Oasis. The total area of blue ice in the study area was estimated to be 106.875 km², approximately 61% of the study area. The WV-2 derived BIR map area presented in this study locally refined the existing BIR map derived using Landsat Enhanced Thematic Mapper Plus (ETM+) and the Moderate Resolution Imaging Spectroradiometer (MODIS)-based mosaic of Antarctica (MOA) dataset by ~31% (~33.40 km²). Finally, we discussed the practical challenges and future directions in mapping BIRs across Antarctica.

1. Introduction

Blue ice regions (BIRs) are areas on ice sheets with negative surface mass balance, where sublimation constitutes the major ablation process and surface albedo is relatively low [1]. They account for around 1% of Antarctica’s surface area [2] and are widely scattered over the continent, although they are generally located in coastal or mountainous regions. In Antarctica, exceptionally dry and windy meteorological conditions favour the formation of large areas of net ablation on the ice sheet, leading to the formation of BIRs [3]. Negative surface mass balance causes fresh surface snow and old firn to be removed primarily through the process of sublimation, resulting in the exposure of dense older ice beneath [1]. Albedo values of various BIRs in different geo-locations around Antarctica have been reported to be in the range of 0.55–0.66 [4]. The lower albedo of BIRs could influence the local surface energy balance and microclimate [5].

Antarctic BIRs are classified into two different types depending on the climate process involved in their formation: (1) melt-induced; and (2) wind-induced [6]. Melt-induced BIRs are located on slopes near coastal areas where surface melt occurs, while wind-induced BIRs occur near mountains or on outlet glaciers where snow is removed by consistent wind year-round. Although melt-induced BIRs constitute a small portion of blue ice, they are an important source of water supply to any surrounding terrestrial habitats, as well as to some Antarctic research stations. Melt-induced BIRs are also utilized as runways for wheeled aircraft [7,8,9,10], such as the runway near the Novolazarevskaya research station in Schirmacher Oasis, East Antarctica. Therefore, one of the main reasons to map BIRs is to find suitable locations for field logistics, as BIRs can facilitate the use of wheeled rather than ski-equipped aircraft.

BIRs have characteristics that make them compelling research locations for planetary geologists, glaciologists, meteorologists, and climatologists. The first glaciological investigations on BIRs were conducted in 1949–1952 by the Norwegian-British-Swedish Antarctic Expedition in Dronning Maud Land (DML) [11], which provided the first description of the surface characteristics of BIRs and discussed their genesis by horizontal compressive forces combined with katabatic winds precluding snow accumulation. A detailed review of the general characteristics of blue ice in Antarctica, the formation and age of BIRs, and their glaciological, meteorological, and climatological features was provided by Bintanja [1].

BIRs are of great potential interest for paleoclimate studies [1,8], although very little paleoclimatic data from BIRs have been published to date [8]. Concentrations of meteorites in Antarctica have been consistently found on, or close to, BIRs [12,13]. Ablation, glacial movement and direct infall are the key natural forces leading to the concentration of meteorites in some BIRs [14]. Seasonal and interannual variations in the spatial extent of BIRs are significantly related to seasonal weather changes [15]; therefore, it is essential to remotely identify BIRs and spatially map their locations. The extent of BIRs is sensitive to different climatic conditions although the details of these relationships are unclear, and there is also feedback from climate change on the process of formation of BIRs. There may be a systematic relationship between temporal variations in the extent of BIRs and their surface elevation, especially those close to nunataks and those undergoing expansion [16]. Some studies have suggested that BIRs can influence the regional surface mass balance and climate and, hence, BIRs could be climate-sensitive [17].

Mapping of the areal extent of BIRs relies heavily on satellite remote sensing data. It is primarily based on multispectral data because of the unique spectral profile of blue ice compared to background features such as snow/rock. The first attempt to map BIRs dates to 1976 when the Japanese National Institute of Polar Research (NIPR) delineated blue ice areas in the Yamato Mountains, Antarctica, using Landsat 1 data [18]. Based on their surface properties, BIRs have been studied using satellite imagery and aerial photographs [19]. Most BIRs are relatively large, so they can be easily observed on satellite images, such as from Satellite Pour l’Observation de la Terre (SPOT), Advanced Very High-Resolution Radiometer (AVHRR), and Landsat. The concept of band (ratio) thresholding has proven to be suitable for mapping BIRs [20,21] and it has, therefore, become a commonly used approach for mapping them in multispectral imagery [2,15,22,23]. Preliminary attempts to map the spatial extent of BIRs have utilized differences in the visible and near infrared (NIR) spectral reflectance of blue ice and snow to map them from medium- and coarse-resolution satellite imagery [2]. Based on the different spectral characteristics (optical satellite imagery) and coherence (SAR imagery) of snow, blue ice, and bare rock, BIRs have been extracted from optical images and coherence maps [24]. Numerous researchers have investigated spectral signatures of blue ice and other cryospheric features (e.g., [23]). BIR mapping has also been attempted with radar imagery [25,26] at the local scale.

A variety of remote sensing (RS) datasets have been used to identify and map BIRs, including aerial photographs and satellite images from Landsat multispectral scanner system (MSS) and TM (e.g., [15]), AVHRR [2,15] and synthetic aperture radar (SAR) images (e.g., [27]). McIntyre [28] mapped BIRs in the Lambert glacial basin by employing Landsat MSS images and estimated their area to be 56,000 km². Swithinbank [29] identified 15 BIRs as potential airfield locations over the Lambert Glacier–Amery Ice Shelf (LGAIS) by utilizing aerial photographs and satellite images. Winther et al. [2] used coarse AVHRR data (1.1 km resolution) acquired in January 1987 and delineated the BIRs of the LGAIS, but the accuracy of the extent estimate was low. Yu et al. [30] mapped spatial extent and geographical distribution of BIRs for the entire LGAIS using Landsat ETM+ images with an accuracy of 92–97%. Brown and Scambos [15] used 56 Moderate-Resolution Imaging Spectroradiometer (MODIS) satellite images to analyse seasonal and interannual variations in the extent of BIRs over the period 1974–2002 near Byrd Glacier on the East Antarctic plateau. All these studies suggest that RS data provided an effective and relatively inexpensive means of detecting and mapping BIRs at various locations in Antarctica.

Nolin et al. [31] used Landsat TM data to map BIRs in the vicinity of the Darwin and Byrd Glaciers, East Antarctica. They also studied surface roughness characterization of blue ice near McMurdo research station using multi-angle imaging spectroradiometer (MISR) data to demonstrate that the blue ice exhibits characteristic angular signatures, which may be used in multispectral classification. Some attempts have been made to map BIRs in Antarctica using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Landsat ETM+, and Ice, Cloud, and Land Elevation satellite (ICESat) laser altimetry to study their different topographical features [30]. ETM+ and MODIS data have been used to map the geographical distribution of BIRs, and supplementary methods were used to study glaciology, meteorites, and surface characteristics of BIRs in Antarctica. Scambos et al. [16] presented an algorithm for mapping snow grain size, which they subsequently employed to estimate areas of selected BIRs in Antarctica using 250 m resolution MODIS bands (band 1: red visible light at 620–670 nm; band 2: near-infrared (NIR) at 841–876 nm). Algorithms implemented for the derivation of BIRs using MODIS and Landsat were based on the increased infrared absorption of ice with grain size [21]. These authors used a normalized difference algorithm for deriving BIRs using Landsat-7 ETM+. The normalized difference snow index (NDSI) was implemented across the Landsat sensor series for long-term analysis of blue ice extent changes. Their work suggested that use of the blue band would potentially further improve the algorithm.

Various other RS methods have been used to discriminate BIRs from firn, snow or exposed rock/nunataks using supervised and unsupervised classification [2], image segmentation using texture [30], and threshold discrimination based on a band ratio [15,16]. The band ratio method has been widely used in studies of snow grain size [32]. Existing normalized difference spectral index ratios (SIRs) utilized in the blue ice mapping application are listed in Supplementary Table S1. Recently, Han et al. [33] developed a Normalized Difference Blue Ice index (NDBI) to map BIRs in the McMurdo Dry Valleys, based on the spectral reflectance characteristics of blue ice, snow, rocks and clouds, using MODIS images. The NDBI, which is defined as (Blue − NIR)/(Blue + NIR), relies on the fact that snow and clouds have a high reflectance in visible and NIR bands. Rivera et al. [34] implemented both manual and automated approaches to map BIRs and crevasses in Patriot Hills, West Antarctica, using ASTER images, GPS, and Radar measurements.

Existing satellite-derived BIR products are either mono-temporal for the pan-Antarctic scale, or multi-temporal on a regional scale. Multitemporal satellite images are necessary for studying seasonal and inter-annual variations of BIRs. However, it is practically impossible to obtain a single temporal satellite image to map BIRs over the entire Antarctic continent, given its vast area (14 × 10⁶ km²). To date, only two BIR mapping applications over the entire Antarctic continent have been implemented. In all other studies, BIRs have been mapped only in limited regional areas. The first map dataset depicting the entire spatial extent of BIRs was generated using AVHRR data with 1.01 km spatial resolution [2]. The second Antarctic-wide study mapped the spatial extent of BIRs and geographical distribution at the unprecedented spatial resolution of 125 m by using the band ratio of Landsat-7 ETM+ images and a snow grain-size image of the MODIS-based Mosaic of Antarctica (MOA) dataset [22] (Figure 1). They used the ratio of band 4 (760–900 nm) and band 7 (2080–2350 nm) to obtain the thresholds for mapping BIRs in Antarctica and estimated a total area of 234,549 km² (1.67% of the area of the continent) during the data acquisition period. Very recently, Hu et al. [35] developed a product of blue ice fraction over Antarctica for the period 2000–2021 using spectral mixture analysis (SMA) performed on moderate-resolution imaging spectroradiometer (MODIS) observations in Google Earth Engine. This approach gives the potential to map BIRs across Antarctica annually.

A comprehensive literature review suggests that the methods currently available for mapping BIRs in cryospheric environments suffer from several important drawbacks. First, most of these methods were developed, evaluated, and validated only for specific regions of Antarctica (e.g., LGAIS, Grove Mountains, Darwin and Byrd Glaciers, McMurdo Dry Valleys, etc.), with most not validated for mapping BIRs in other areas. Second, most of them used either medium- or coarse-resolution (SPOT, AVHRR, MODIS, Landsat and ASTER) optical satellite RS data. As the spatial resolution of an image increases, it is possible to detect smaller blue ice patches and to improve mapping accuracy. Due to the coarser pixel size of MODIS, ASTER and Landsat-7 data, BIR mapping is restricted to the larger BIRs. Furthermore, most of the methods attempted have not been thoroughly compared to existing traditional image-processing methods. In this situation, snow patches of varying sizes and shapes which are identified on the BIRs cannot be eliminated using coarse/medium resolution satellite imagery. Therefore, it is necessary to extract BIRs utilizing high-resolution RS data to refine their spatial distribution. Third, to our knowledge, there is no automatic or semi-automatic BIR extraction method that uses VHR satellite RS data for the Antarctic environment. Hence, this is the first comprehensive study applying a semi-automatic BIR mapping approach developed using VHR WV-2 data in the Antarctic environment.

Based on the literature review and the limitations of current geospatial methods discussed above, we focused our study on mapping BIRs in the vicinity of the Schirmacher Oasis, Princess Astrid Coast, Dronning Maud Land, East Antarctica, using a normalized difference algorithm and pixel-based supervised machine learning methods. We performed semi-automatic extraction of BIRs by employing four types of feature extraction approach: (a) a customized NDBI approach; (b) a spectral processing (SP) or matching approach; (c) pixel-wise supervised classifier (PSC); and (d) a target detection (TD) approach. The main goal of this study was to evaluate the potential of pixel-based band ratio and supervised machine learning methods for an effective mapping of BIRs using WV-2 data. Our study focused on the following objectives: (a) designing a customized NDBI approach to extract BIRs; (b) comparing the performance of supervised feature extraction algorithms with the newly developed customized NDBI approach using visual analysis and statistics; and (c) assessing the distinctive 8-band acquisition capability of WV-2 data by using an array of semi-automatic mapping methods to compare their reliability in BIR extraction. The precise, high-resolution, and refined BIR map derived from this research will provide the baseline information required for future climate change analysis of the BIRs in this region.

2. Study Area and Geospatial Data

2.1. Extent of the Study Area

The Schirmacher Oasis (SO) is a ∼35 km long and up to ∼3 km wide ice-free plateau (Figure 1) between 70°43′50″ and 46′40″S, and 11°22′40″ and 54′20″E, with more than 100 freshwater lakes. The SO is located on the Princess Astrid Coast of Dronning Maud Land, East Antarctica. It has a low-lying hilly topography, with a maximum altitude of 228 m above sea level (asl) and an average altitude of ~100 m asl. The SO is located between the edge of the East Antarctic Ice Sheet (EAIS) and the Novolazarevskaya Nivl iceshelf. The plateau of the SO provides a barrier to the northwards running ice stream. South of the SO, the EAIS reaches a height of 1500 m asl at a distance of 50 km from the oasis. Several nunataks are exposed above the ice sheet in this slope area.

The general climate of the SO is relatively mild compared with the overall Antarctic climate. Air temperature ranges from −7.7 to +8.2 °C during mid-summer (December–February) resulting in abundant meltwater. January is the warmest month (monthly mean air temperature +0.7 °C, maximum +8.2 °C) while August is the coldest (monthly mean air temperature −16.3 °C, minimum −35.5 °C), with an average wind velocity of about 9.7 ms⁻¹ and 264.5 mm annual precipitation, mostly as snow [36]. Atmospheric pressure, air temperatures and air humidity are greatest in summer, with comparatively weak winds, rare snowfall, and insignificant total precipitation. Typically, rapid melting of snow/ice and drainage of meltwater from the SO to the iceshelf is observed. In general, katabatic winds over the EAIS increase towards the coast and become more directed down the fall line of the topography. This induces a divergence of drifting snow transport that serves as a negative term in the surface mass balance (erosion). The blue ice extent has been relatively stable near exposed nunataks in Dronning Maud Land [37].

2.2. Spatial Data Used in This Study

A list of satellite and ground reference datasets and their usage in present study is given in

Table 1. Satellite and ground reference datasets and their specific usage in this study. A variety of spatial datasets (both satellite data and ground reference data) were utilized at different stages of the study. The potential usage of these datasets includes: (1) SE (semi-automatic feature extraction): The datasets are used as input for semi-automatic extraction methods; (2) VI (visual interpretation): The datasets are used for visual interpretation to support digitization, to interpret BIRs, and to ensure spatial snow/ice cover changes are captured; (3) MD (manual digitization): The datasets are used in an ArcGIS platform to carry out digitization; (4) SD (supplementing digitization): The datasets are used to cross-check and verify the manual digitization against various ground reference and satellite images; and (5) DEA (digitization error analysis). The datasets are used to test positional errors in manual digitization. The notation ✓ (🗴) denotes used (not used) in this study.

Dataset	Source	Temporal Range (dd/mm/yy)	Utilization of Datasets in the Present Study
			SE	VI	MD	SD	DEA
Worldview-2 MSI and PAN	DigitalGlobe	5 February 2012	🗴	✓	🗴	✓	🗴
Worldview-2 PAN-sharpened image (0.5 m)	Processed	5 February 2012	✓	✓	✓	✓	✓
Google Earth Images (GE)	Google	31 December 1999 to 26 November 2013	🗴	✓	🗴	✓	🗴
BIR map derived from MODIS and Landsat ETM+	Hui et al. [22]	1999–2004	🗴	✓	🗴	✓	🗴
DGPS Surveying	InSEA	2008–2015 (September–March)	🗴	✓	✓	✓	✓

.

2.2.1. Usage of Satellite Data

We employed radiometrically corrected, georeferenced, ortho-rectified 16-bit standard level 2 (LV2A) WV-2 multi-sequence data, acquired on 5 February 2012, at an off-nadir angle of 17.02° over the SO. Based on the results of MODIS-based BIR mapping [15], we selected the image as close to mid-summer as possible. The data consists of four tiles of 8-band multispectral (MS) data acquired by visible-infrared (V-NIR) WorldView-110 camera and a panchromatic (PAN) image. These tiles were spatially mosaicked to generate a single continuous image. The projection and datum of the SO images were geo-registered with UTM zone 32S and WGS 1984, respectively. The WV-2 image covers an area of ~172 km² including various land-cover features (snow, ice, blue ice, rocks, lakes, permafrost, vegetation, etc.), and terrain consisting of flat areas up to mountains of 600 m altitude. Since the satellite image was acquired at the height of the austral summer, solar radiation could melt surface snow and ice on clear days. Every spectral band is closely concentrated on a specific range of EMR wavelengths sensitive to a specific target feature on the ground. A WV-2 MS image (MSI) consists of four traditional bands: Band 2, Blue (450–510 nm), Band 3, Green (510–580 nm), Band 5, Red (630–690 nm), and Band 7, NIR-1 (770–895 nm), and four new bands: Band 1, Coastal (400–450 nm), Band 4, Yellow (585–625 nm), Band 6, Red edge (705–745 nm), and Band 8, NIR-2 (860–1040 nm).

2.2.2. Ground Truthing and Supplementary Data

The ground reference datasets utilized to support the semi-automatic extraction of BIRs were retrieved from the Indian scientific expeditions to Antarctica (InSEA), a continent-wide BIR map derived from MODIS and Landsat-7 ETM+ [22], historical GE images (Supplementary Figure S1), and PAN-sharpened WV-2 images. Historical GE images acquired during the austral summer were used to cross-check the spatiotemporal variation of blue ice extent in the study area.

Extensive field surveys have already been carried out by teams participating in the InSEA since 2005. A database consisting of BIRs in the SO region was generated by manual digitization. PAN-sharpened images (0.5 m) were visualized in ArcGIS 10 at several scales using various band combinations of WV-2 data: 7-4-2, 8-7-2, 6-3-2, 5-3-2 and 7-3-2. However, we could easily recognize the pattern, texture, shape, and size (dimension) of BIRs and discriminate them against background features using a 5 (red)-3 (green)-2 (blue) band combination and 1:500 scale, so this combination was used for digitization of BIRs. Each BIR polygon was then edited to eliminate digitization errors. Quality control involved cross-checking the BIR polygon boundaries to ensure accurate and consistent interpretation and digitization. Then, the positional accuracy of the digitized BIR database was calculated using an independent source of 20 accurate differential global positioning system (DGPS) point locations [38]. The RMSE of digitization yielded ~30 cm (less than 1 pixel) when compared with the DGPS source. The final manually digitized reference map depicting distribution of BIRs is shown in Figure 2. The manually digitized 14 BIR tiles were extensively ground-surveyed, confirmed and cross-checked using DGPS ground reference data obtained from field campaigns, the MODIS and ETM+ derived BIR map and multitemporal GE images.

A field survey of 14 tiles was carried out in the Real-time kinematic (RTK) mode using a Leica Viva DGPS surveying unit. Since these field datasets were collected from 2011 to 2015 (September–February), they were processed to conform to the 2012 (February) images. Visual examination was conducted to eliminate those BIR patches that were visible in the existing WV-2 imagery because they may have been covered by seasonal snow or supraglacial debris. Therefore, only 12 tiles (Figure 1) mapped using the WV-2 image acquired on 5 February 2012 were considered in the present analysis to reduce the errors associated with temporal changes in BIRs. Those 12 tiles mapped using DGPS recorded during January-February were overlaid on the WV-2 PAN-sharpened image and checked against the digitized reference database. After a careful visual analysis, very few (~8%) manually digitized BIRs that had higher boundary variations (typically > 1 pixel, i.e., ~50 cm) compared with the ground-surveyed BIR boundaries were corrected manually.

A GIS–compatible shapefile of the BIRs was generated using ArcGIS 10. The surface area of BIRs and Non-BIRs (snow/supraglacial debris/streams) for all 12 tiles was calculated using GIS routines. A list of the 12 tiles and their reference areas is given in

Table 2. Sample data summarizing the reference BIRs (manually digitized and cross-verified with ground surveying) for the 12 tiles and their respective non-BIR areas (snow/ice/supraglacial debris, streams etc.). Topographic distribution of tiles is highlighted with coloured text and background; (i) tiles on comparatively lower elevation areas are highlighted in violet text (Tiles 8, 9, 12), (ii) tiles on comparatively higher elevation areas are highlighted with green text (Tiles 2, 4, 6), (iii) tiles on intermediate elevation areas are highlighted with red text (Tiles 1, 3, 5, 7, 10, 11), (iv) tiles including shadow-prone areas are highlighted with underlined text (Tiles 8, 9), and (v) Tiles including non-shadowed areas are highlighted with blue background (Tiles 1, 2, 3, 4, 5, 6, 7, 10, 11, 12). Spatial distribution of tiles based on surface area of BIRs is highlighted with bold italics underlined (<65% blue ice, comparatively low blue ice, Tiles 6, 8, 9), bold underlined (>80% blue ice, relatively high blue ice, Tiles 1, 3, 4, 5, 7), and bold italics (65–80% blue ice, Tiles 2, 10, 11, 12).

Tile No.	Total Area (m2)	Reference BIRs		Reference Non-BIRs
		(m2)	%	(m2)	%
1	202,392.18	189,504.78	93.63	12,887.40	06.37
2	202,119.31	133,320.08	65.96	68,799.23	34.04
3	202,257.09	164,433.85	81.30	37,823.24	18.70
4	201,547.57	170,673.36	84.68	30,874.21	15.32
5	201,630.23	169,468.74	84.05	32,161.49	15.95
6	201,608.05	127,106.45	63.05	74,501.60	36.95
7	201,750.89	190,304.24	94.33	11,446.65	05.67
8	202,171.24	104,841.94	51.86	97,329.30	48.14
9	201,947.31	104,285.23	51.64	97,662.08	48.36
10	201,588.89	146,478.47	72.66	55,110.42	27.34
11	201,444.93	137,480.92	68.25	63,964.01	31.75
12	201,520.04	137,752.88	68.36	63,767.16	31.64

. Sample reference data for Tiles 1 and 2 are depicted in Figure 2. In the present analysis, the 12 tiles were selected on the basis of the following criteria (Figure 1 and Supplementary Figure S2): (a) geographical regions represented by the selected 12 tiles were accessible for conducting field surveys; (b) to ensure inclusion of BIRs of varying spatial extent across these tiles to support the robustness of the analysis; (c) the 12 tiles were spatially well-distributed to ensure unbiased accuracy analysis; (d) the 12 tiles were selected to include various elevations to ensure the robustness of the analysis against topographical errors; and (e) varying spatial extents of the confounding spectral signatures from non-target areas (non-BIRs) were included over the 12 tiles to test the robustness of the analysis against various non-BIR features.

3. Methods and Analyses

The data processing protocol is shown in Figure 3. The steps consisted of three blocks as described below: (a) data pre-processing; (b) blue ice mapping using feature extraction (FE); and (c) quantitative evaluation of accuracy.

3.1. Pre-Processing of Satellite Data

Two procedures were implemented for data pre-processing: (a) WV-2 satellite data calibration; and (b) atmospheric corrections. The conversion of DN to at-sensor radiance and subsequently to top-of-atmosphere (TOA) reflectance images was conducted using spectral transmission of the telescope, the spectral quantum efficiency of the detectors, and other calibration parameters provided in the metadata file by using ENVI 5.x. The top-of-atmosphere (TOA) reflectance images were converted to at-surface reflectance using ATmospheric and Topographic CORrection version 3 (ATCOR-3 algorithm) (Supplementary Figure S3) [39] coupled with a digital elevation model (DEM) constructed using synergetic merging of ICESat-1 and DGPS data for the SO [40]. Parameters used for ATCOR-3 processing are summarized in Supplementary Table S2. Finally, to create a sharpened image at 0.50 m resolution, the multiband MS image was PAN-sharpened using the Gram-Schmidt (GS) method, which has been tested and proved to be comparatively superior and consistent for sharpening WV-2 data [41].

3.2. Blue Ice Detection and Mapping Using Feature Extraction

This step was conducted using four semi-automatic FE approaches: (a) a customized NDBI approach (cNDBI); (b) a spectral processing or matching (SP) approach; (c) a target detection (TD) approach; and (d) pixel-wise supervised classification (PSC). All satellite data were processed using a procedure developed in ENVI software (L3Harris Geospatial, Colorado, Boulder) and Interactive Data Language (IDL; L3Harris Geospatial, Colorado, Boulder).

3.2.1. Customized NDBI Approach

The physical basis for most of the RS methods developed for mapping BIRs can be summarized as follows. (1) A lower surface albedo of blue ice (0.5–0.6) compared to snow (0.8–0.9) causes BIRs to exhibit a darker surface in contrast to the surrounding snow [1,42]; (2) blue ice has minimum absorption at 470 nm. BIRs can also be mapped by comparing visible bands with a near- or middle-infrared channel [16]. The band ratio method is principally based on the characteristic differences in spectral reflectance rendered by the blue ice, which can be used for distinguishing blue ice from non-target classes (snow, frozen lakes, supraglacial lakes, moraines, or exposed rock). As shown by its spectral curve (Supplementary Figure S4), blue ice strongly absorbs solar radiation in the red and infrared wavelengths and reflects in the blue portion of the spectrum, which explains its bluish appearance. In comparison with blue ice, fresh snow has much higher reflectivity (and hence reduced solar radiation absorption), particularly at red and near-infrared (NIR) wavelengths of the spectrum permitting easier discrimination of blue ice from snow [23]. A comparison of spectral curves of blue ice and snow indicates that the red and NIR bands have the optimal spectral separability to discriminate them. Due to their apparent spectral differences (strong absorption in the visible and NIR), blue ice and exposed rocks can be visually separated from snow background on the colour composite of WV-2 bands, particularly on the false colour composite of green, red and NIR bands.

The existing band ratios make use of the shortwave infrared (SWIR), red, blue and NIR bands to extract BIRs [16,33]. However, for WV-2 data, the blue ice indices (e.g., NDBIs) were revisited to incorporate the new bands and design customized NDBIs for blue ice mapping and to provide a wider context to the analysis. Hence, in addition to existing band ratios, we proposed various customized NDBIs in this study to fully utilize the add-on WV-2 bands. Our unique set of customized NDBIs was designed based on the minimum redundancy maximum relevance (mRMR)-based band scoring/ranking (Supplementary Figure S5) and visual analysis (Supplementary Figure S4) of the spectral responses of various BIRs present on the WV-2 image. It is a time consuming-process for WV-2 data, because of its high spatial and spectral resolution compared with other datasets.

In this study, the new spectral bands (Coastal, Yellow, Red Edge, NIR1, and NIR2) and traditional bands (Blue, Green, Red) were ranked for maximum and minimum response values relevant to a classification; bands were also ranked using the mRMR criterion. We selected target blue ice regions in MSI and then found the combination of spectral bands that differentiated the target class. Spectral bands were ranked based on their ability to represent the target class (mRMR criterion), which was visually cross-verified. The bands ranked according to mRMR were empirically evaluated by observing the spectral profiles for BIRs. We evaluated nearly 1000 spectral profiles from the blue ice class to ensure the maximum and minimum response band. These bands were normalized to yield the normalized difference target class index (NDCI), in this case, NDBI. The detailed procedure of customizing NDBIs is described in Figure 4. The spectral profile analysis and mRMR-based scoring depicted in Supplementary Figure S4 show that the reflectance values range from 40–100% for blue ice and 95–100% for snow.

The reflectance values for blue, coastal, green, and yellow bands are comparable for most of the profiles. The maximum reflectance rendered by these four bands can be used to identify the blue ice on the image. The customized NDBIs proposed for establishing effective blue ice mapping methods using WV-2 images are listed in

Table 3. Customized NDBI devised as an effective semi-automatic extraction tool for blue ice mapping application. A list of threshold ranges defined for the customized NDBIs to effectively extract blue ice areas is also given.

NDBI	NDBI Model	Mathematical Expression	Threshold Range
NDBI-1	NDBI(3−7/3+7)	$\frac{\mathrm{Green}-\mathrm{NIR}1}{\mathrm{Green}+\mathrm{NIR}1}$	0.83−0.95
NDBI-2	NDBI(3−8/3+8)	$\frac{\mathrm{Green}-\mathrm{NIR}2}{\mathrm{Green}+\mathrm{NIR}2}$	0.87−0.92
NDBI-3	NDBI(4−7/4+7)	$\frac{\mathrm{Yellow}-\mathrm{NIR}1}{\mathrm{Yellow}+\mathrm{NIR}1}$	0.84−0.93
NDBI-4	NDBI(4−8/4+8)	$\frac{\mathrm{Yellow}-\mathrm{NIR}2}{\mathrm{Yellow}+\mathrm{NIR}2}$	0.85−0.96

. Sample spectral profiles (Supplementary Figure S4) of blue ice reveal that coastal, blue, yellow, and green bands have maximum mRMR score for blue ice and can be explored for customizing NDBI. Following this criterion, NDBIs were designed using the relative responses of blue ice in the spectral space (Supplementary Figure S4). The NDBI is used to identify areas where BIRs are the dominant background or foreground material. Based on the spectral signature analysis, SWIR and NIR spectral channels can be most effectively used to characterize the difference in reflection for BIRs. The existing shortwave infrared (SWIR) band-based NDSI is analogous to the traditional NDVI (NIR-Red normalization) [43] which is useful for discriminating and mapping snow/ice from cumulus clouds. The NDSI is usually used for mapping snow cover [44], based on the high and low reflectance of snow in the visible (Green) and SWIR regions, respectively, yielding the ratio (Green − SWIR)/(Green + SWIR). Brown and Scambos [15] implemented a NDSI based on NIR and green bands of Landsat-7 ETM+ to map blue ice. Since there is no SWIR band available on WV-2, it was necessary to customize the traditional NDSI to generate a blue ice index for this study. However, there are characteristic differences in the response values for BIRs in the Yellow–NIR-1, Yellow–NIR-2, Green–NIR-1 and Green–NIR-2 bands, which have been exploited to identify blue ice (Table 3). This modification in traditional NDSI and NDVI increases the sensitivity of the NDBI algorithm relative to NDVI because of the higher, more uniform reflectivity for snow and ice in the green and yellow channels. Using the blue or coastal band would potentially improve the algorithm still further, but scattered light from the atmosphere is a significant problem in the blue channel; therefore, green and yellow bands were selected to design the NDBI. We note that the NDBI does not maximize the response of blue ice, but best discriminates the target blue ice class from other regions in the image.

Drawing from a sample of the blue ice within the scene, a consistent and unique difference between the customized NDBI bands was observed, which was constant for all BIRs (Supplementary Figure S4). All detections were validated by cross-referencing the signatures in corresponding regions in the true and false colour composites with the detections in the NDBI. This is a novel method for determining BIRs without using a SWIR band. NDBI produces a single grayscale image, where blue ice appears bright. Some other features, such as frozen lakes, epishelf lakes and meltwater on the surface of ice also appear bright. The possibility of supraglacial lakes having a spectral response like that of blue ice in the study area is remote because there are no supraglacial lakes on the WV-2 imagery of the study region that have an analogous spectral response to blue ice. Moreover, the satellite data were captured with minimum cloud cover, and blue ice was the only target that appeared bright on the NDBI image.

Thereafter, the customized NDBIs were calculated against the values in the reflectance cube, and an output image was generated for each customized NDBI. Because of increased solar radiation meltwater accumulated on the surface of white ice, on the surface of frozen or semi-frozen lakes and on epishelf lakes, which produced a spectral response similar to NDBI. To discriminate between blue ice and meltwater on the surface of snow/lakes, it was essential to threshold the NDBI value so that the entire region was classified as the target class (blue ice) or a non-target class (snow, white ice, meltwater, frozen lake, epishelf lake). Since the BIRs exhibited NDBI values ranging from 0.84 to 0.96, a scene-dependent threshold was defined and used to discriminate between blue ice and non-blue ice pixels. Threshold values used in this study were empirically evaluated by manually scrutinizing the most obvious 1000 blue ice pixels from the NDBI images obtained from the PAN-sharpened images. Analyses were performed at the same pixel locations for all NDBI images derived from PAN-sharpened images to verify unbiased accuracy analysis. Table 3 summarizes the thresholds defined for the customized NDBI methods. A detailed description of the thresholding process is given in Supplementary Material S1. Pixels with higher NDBI values than the local threshold were coded as 1 (target blue ice class pixels), while pixels with a lower NDBI value were coded as 0 (non-target class pixels). After classifying the image based on each customized NDBI, semi-automatically extracted BIRs were vectorized to compute the area. Customized NDBI-based extracted blue ice area (m²) and the bias (m²) with respect to the reference digitized blue ice area (m²) for all 12 tiles are summarized in

Table 4. Sample data summarizing the reference blue ice area (manually digitized and cross-verified with ground survey) for the 12 tiles containing blue ice and their respective blue ice areas extracted by the customized NDBI approach using GS-sharpened images. Negative bias values representing overall overestimation of blue ice are highlighted with bold italics and underlined text. Outliers are highlighted with *.

	NDBI-1		NDBI-2		NDBI-3		NDBI-4
Tile No.	Extracted Area (m2)	Bias in Area (m2)	Extracted Area (m2)	Bias in Area (m2)	Extracted Area (m2)	Bias in Area (m2)	Extracted Area (m2)	Bias in Area (m2)
1	189,540.73	−35.95	189,741.59	−236.81	189,913.57	−408.79	189,969.43	−464.65
2	132,971.17	348.91	133,071.89	248.19	133,193.03	127.05	133,293.27	26.81
3	164,849.04	−415.19	164,901.71	−467.86	164,831.11	−397.26	164,801.11	−367.26
4	170,279.86	393.50	170,079.25	594.11	170,411.25	262.11	170,411.25	262.11
5	169,810.47	−341.73	169,972.57	−503.83	170,802.57	−1333.83	170,831.73	−1362.99
6	126,907.83	198.62	126,986.98	119.47	127,031.04	75.41	126,931.04	175.41
7	189,928.71	375.53	190,011.38	292.86	189,901.83	402.41	190,621.57	−317.33
8	103,984.90	857.04	103,289.61	1552.33	103,937.55	904.39	103,326.68	1515.26 *
9	103,257.66	1027.57	103,097.47	1187.76	103,203.41	1081.82	103,089.41	1195.82
10	146,967.39	−488.92	147,167.39	−688.92	145,716.82	761.65	146,599.38	−120.91
11	137,598.07	−117.15	137,971.44	−490.52	137,991.03	−510.11	138,381.74	−900.82
12	137,421.31	331.57	137,861.83	−108.95	138,287.58	−534.70	137,399.43	353.45
|Avg.|	147,793.10	410.97	147,846.09	540.97	147,935.07	566.63	147,971.34	588.57
RMSE		491.65		682.58		675.21		768.46

.

3.2.2. Spectral Processing (SP) or Matching–Based Extraction Approach

Spectral processing/matching methods extract the target features from MS imagery based on the target feature’s spectral characteristics. The supervised workflows were used with an initial approximation in terms of Regions of Interest (ROIs) as reference spectra that include representative pixels of BIRs. Spectral processing methods determine the inherent spectral similarity between input spectral bands and reference spectra to generate an output product showing pixels with similar spectral properties clumped into target and non-target classes. A set of four SP methods were conducted using Spectral Processing Exploitation and Analysis Resource (SPEAR) workflow tools and mapping methods (ENVI 5.x), which streamline methods [matched filtering (MF), mixture tuned matched filtering (MTMF), mixed filtering/spectral angle mapper (MF/SAM) and principal component analysis (PCA)] for mapping BIRs using WV-2 PAN-sharpened data. BIRs covered with varying extents of snow showed a wide range of shades of blue colouration, ranging from light blue to dark blue. The major influences came from topography of the study region and varying amounts of snow sporadically scattered in many regions. Because of this difference in spectral information, unsupervised information extraction methods were not used for automated BIR mapping. Since one of our objectives was to extract BIRs semi-automatically, we processed the imagery using supervised workflows based on knowledge of BIR locations within the landscape.

3.2.3. Target Detection (TD) Approach

TD methods work on the principle of extracting target objects based on spectral characteristics of reference training spectra (seed points) of target features and suppressing background noise using spectra of non-target features. These methods provide a rapid means to maximize the response of the known target spectra and suppress the response of the composite unknown background targets. The initial approximations or spectral signatures for supervised workflows are provided in terms of target ROIs and non-target ROIs. Target detection tools (ENVI 5.x) were used to perform supervised image processing tasks into workflows [constrained energy minimization (CEM), adaptive coherence estimator (ACE), orthogonal subspace projection (OSP), and mixed tuned target-constrained interference-minimized filter (MT-TCIMF)] to extract blue ice features. We used the TD workflow that couples a minimum noise fraction (MNF) transformation for improved TD for multispectral data [45]. Additionally, for MT-TCIMF, the MNF-transformed image fuses the mixture-tuned method with the TCIMF method. Hence, the MNF was used to resolve the intrinsic spectral dimensionality of the WV-2 image to support subsequent processing. The MNF transformation served two purposes: (a) MNF was used to extract target features and/or reduce dimensions for the generic TD methods, such as ACE, CEM, and OSP; and (b) for the mixture-tuned method (MT-TCIMF), applying the MNF transformation resulted in the isotropic variance noise, which was used to compute infeasibility value. By adding an additional infeasibility band, mixture-tuned techniques improved the detection results by reducing the number of false spectral signals.

3.2.4. Pixel-Wise Supervised Classification Approach

We used four popular and distinctly different pixel-based classification methods to classify the WV-2 PAN-sharpened data: maximum likelihood classifier (MXL), support vector machine (SVM); neural network classification (NNC); and SAM in ENVI 5.2. The robustness of these classifiers in mapping land-cover features using VHR image is discussed elsewhere [46]. The Gram Schmidt (GS)-sharpened WV-2 image was classified into target (BIRs) and non-target (meltwater, snow, landmass, lakes) areas. We limited our analysis to these two classes (BIRs and non-BIRs) to reduce the time needed to produce the final map of BIRs.

3.2.5. Practical Execution of Image Processing Routines and Post Processing Corrections

The FE methods described in this section are based on different underlying principles. To compare these methods objectively, we kept the input ROIs (training samples) constant for all methods. For the sharpened image, we selected 1000 regions of interest (ROIs) for non-target and target (BIRs) classes as training pixels using spectral signature observations. The training datasets were randomly selected based on manual interpretation with the GS-sharpened image, field checking and using multi-date VHR historical images from GE. Individual reference spectra (end members) for the blue ice class were collected as ROIs over spatially and spectrally homogeneous targets on WV-2 PAN-sharpened images. End member or ROI determination was done by carefully considering the blue ice unit size to be identified and sampled on the image. By selecting such separable samples with the ROI tool, we reduced the problem of mixed pixels. Spectral signatures or reference spectra for each training class were defined by calculating image statistics of ROIs. Finally, FE methods were executed using the common input ROIs (training samples) for all three approaches, i.e., SP, TD and PSC.

BIRs were identified using 16 pixel-based methods and the resulting classified binary images (0: non-BIRs, 1: BIRs). To correct possible classification errors and to suppress or reduce the random noise, a number of post-classification operations were performed, including majority filtering and mathematical morphology operations (dilation, erosion, fill, trim) [47]. To reduce a ‘salt-and-pepper’ effect and produce more realistic maps, a median filter with 3 × 3-pixel windows was applied. This filter is less sensitive to extreme changes in pixel values as compared to linear techniques; it can remove salt-and-pepper noise and preserve useful details in images without significantly reducing the sharpness of an image. Finally, after classifying the image into the target class, i.e., blue ice, using the SP, PSC or TD approaches, the semi-automatically extracted BIRs were vectorized to calculate the blue-ice area derived from the individual method. Thereafter, we compared these calculated areas with manually digitized reference blue ice areas, and we evaluated the statistical significance based on the accuracy assessment.

3.3. Accuracy Assessment

The results of the extraction methods were assessed using visual and statistical analyses. We simultaneously overlaid blue ice portions of the original PAN-sharpened image, the resulting semi-automated extracted BIRs, and the manually digitized reference BIR shapefile, to perform identical operations on these images within a single image window. We randomly picked a blue-ice polygon feature from the image to check the accuracy of the semi-automatic extraction method by observing the variation in compactness and completeness of the outline of the extracted blue ice polygon compared with the outline of the manually digitized reference blue ice polygon. Statistical analyses were conducted in three different steps: (i) comparing tile-wise accuracy analysis for 12 manually digitized and ground surveyed tiles; (ii) comparing areas of extracted BIRs against manually digitized BIRs for the entire study region; and (iii) comparing the extracted BIRs against the existing BIR map. Root mean square error (RMSE) and bias values between manually digitized (reference) and semi-automatically extracted blue ice extent were estimated. Outliers in bias and RMSE were calculated based on interquartile range (IQR) (i.e., Tukey’s fences) in which Q₁ and Q₃ are the lower and upper quartiles respectively, and outlier was defined as any value outside the range: [Q₁ − k (Q₃ − Q₁), Q₃ + k(Q₃ − Q₁)], where k = 1.5. To test the robustness of the customized NDBI approach in extracting BIRs from topographically varied areas, we organized the 12 tiles in three groups: (i) tiles on low elevation areas (tiles 8, 9 12) (LET); (ii) tiles on high elevation areas (HET) (tiles 2, 4, 6); and (iii) tiles on moderate elevation areas (MET) (tiles 1, 3, 5, 7, 10, 11). Statistical parameters were calculated separately for these three groups of tiles to compare the results of blue ice extraction.

4. Results

The performance of the customized NDBI-based semi-automatic extraction approach was compared with the supervised semi-automatic FE approaches by computing RMSE and bias of the extracted blue ice areas. The performance of the 16 feature extraction algorithms (grouped into four approaches) was also compared for their applicability to extract BIRs.

4.1. Comparing Tile-Wise Accuracy for 12 Manually-Digitized and Ground-Surveyed Tiles

A pre-digitized database of BIRs of the 12 extracted tiles in the form of shapefiles was utilized as a reference for assessing the accuracy of the semi-automatic extraction methods (Table 2). The vectorized area of BIRS from the automated extraction procedures was compared with BIRs in the manually digitized database. For quantifying the uncertainty of our analysis, we used the RMSE because it is a consistent and statistically significant indicator of accuracy. The RMSE statistics were calculated to evaluate the accuracy of extraction and misclassification based on false signals. As the resolution of the WV-2 PAN-sharpened image is 0.5 m, each misclassified pixel or false signal can introduce a RMSE of 0.25 m² in the semi-automatically extracted blue ice area. Mathematical expressions for calculating various parameters of accuracy analyses are shown in Supplementary Material S1.

4.1.1. Approaches for Semi-Automatic Extraction of BIRs

The average and total RMSE values for each FE approach (column) were calculated to judge the overall stability and consistency of each approach (

Table 5. Quantitative evaluation (bias, % and m²) of various feature extraction methods for extracting blue ice from the 12 tiles analysed. The lowest (italics) and highest (bold) values in each local column are highlighted. The column-wise total and average values are bolded and underlined. Total bias represents the total misclassified blue ice area, positive bias represents the underestimated area of the total misclassified blue ice area and negative bias represents the overestimated area of the total misclassified blue ice area. Local outliers (w.r.t. four methods of the respective feature extraction approach) are highlighted with *. RMSEt refers to the error values pertaining to the 12 tiles while RMSEm refers to error values pertaining to the four methods in each respective approach.

		Total Bias (m2)	Total Bias (%)	Positive Bias (m2)	Positive Bias (%)	Negative Bias (m2)	Negative Bias (%)	RMSEt (m2)
NDBI-1	Present work	4931.68 *	0.28 *	3532.74	0.20	−1398.94	0.08	491.65 *
NDBI-2	Present work	6491.61	0.37	3994.72 *	0.22 *	−2496.89	0.14	682.58
NDBI-3	Present work	6799.53	0.38	3614.84	0.20	−3184.69	0.18	675.21
NDBI-4	Present work	7062.82	0.40	3528.86	0.20	−3533.96	0.20	768.46
	Average	6321.41	0.36	3667.79	0.21	−2653.62	0.15	654.48
	RMSEm	6375.33	0.36	3672.80	0.21	2775.91	0.16	662.21
MT-TCIMF	[48]	9286.29	0.52	4423.17	0.25	−4863.12	0.27	860.71
OSP	[49]	11,063.73	0.62	5462.31	0.31	−5601.42	0.32	1035.55
ACE	[50]	10,167.08	0.57	5872.27	0.33	−4294.81	0.24	941.54
CEM	[51]	12,551.57	0.71	6837.42	0.39	−5714.15	0.32	1113.45
	Average	10,767.17	0.61	5648.79	0.32	−5118.38	0.29	987.81
	RMSEm	10,834.58	0.61	5714.79	0.32	5150.79	0.29	992.40
MTMF	[52]	12,556.91	0.71	6064.41	0.34	−6492.50	0.37	1119.16
MF/SAM	[53]	12,899.18	0.73	6822.03	0.38	−6077.15	0.34	1152.36
MF	[54]	15,410.72	0.87	8701.81	0.49	−6708.91	0.38	1361.74
PCA	[55]	19,144.35	1.08	10,500.32	0.59	−8644.03 *	0.49 *	1675.11
	Average	15,002.79	0.84	8022.14	0.45	−6980.65	0.39	1327.09
	RMSEm	15,232.05	0.86	8205.11	0.46	7050.05	0.40	1345.44
SVM	[56]	22,520.21	1.27	12,598.20	0.71	−9922.01	0.56	1976.81
SAM	[57]	25,156.87	1.42	12,910.34	0.73	−12,246.53	0.69	2179.90
NNC	[46]	27,370.25	1.54	17,488.27 *	0.98 *	−9881.98	0.56	2503.42
MXL	[58]	26,526.68	1.49	11,605.36	0.65	−14,921.32	0.84	2377.57
	Average	25,393.50	1.43	13,650.54	0.77	−11,742.96	0.66	2259.43
	RMSEm	25,459.89	1.43	13,837.59	0.78	11,923.96	0.67	2268.24
	Total Avg.	14,371.22	0.81	7747.32	0.44	−6623.90	0.37	1307.20
	Total RMSEm	16,110.98	0.91	8731.42	0.49	7518.65	0.42	1447.28

). In addition, RMSE (m²), bias (% and m²), positive bias or total underestimated blue ice area (% and m²), negative bias or total overestimated blue-ice area (% and m²) values for each extraction method were calculated to assess the performance of each method within a respective extraction approach (Table 5 and

Table 6. Quantitative evaluation (misclassified number of pixels, % and #) of various feature extraction methods for extracting blue ice from the 12 analysed tiles. The lowest (italics) and highest (bold) values in each local column are highlighted. The column-wise total and average values are bolded and underlined. [MP: Misclassified Pixels, OP: Overestimated pixels, UP: Underestimated pixels]. Global outliers (w.r.t. the 16 methods) are highlighted with (′) and local outliers (w.r.t. the four methods in each respective approach) are highlighted with *.

		MP (#)	OP (#)	UP (#)	OP (%)	UP (%)
NDBI-1	Present work	19,727 *	14,131	5596	71.63 ′	28.37 ′
NDBI-2	Present work	25,966	15,979 *	9988	61.54	38.47
NDBI-3	Present work	27,198	14,459	12,739	53.16	46.84
NDBI-4	Present work	28,251	14,115	14,136	49.96	50.04
	Average	25,286	14,671	10,615	59.07	40.93
	RMSEm	25,501.17	14,691.06	11,103.87	59.67	41.78
MT-TCIMF	[48]	37,145	17,693	19,452	47.63	52.37
OSP	[49]	44,255	21,849	22,406	49.37	50.63
ACE	[50]	40,668	23,489	17,179	57.76	42.24 *
CEM	[51]	50,206	27,350	22,857	54.48	45.53
	Average	43,069	22,595	20,474	52.31	47.69
	RMSEm	43,338.16	22,859.23	20,603.21	52.46	47.86
MTMF	[52]	50,228	24,258	25,970	48.3	51.70
MF/SAM	[53]	51,597	27,288	24,309	52.89	47.11
MF	[54]	61,643 *	34,807	26,836	56.47	43.53
PCA	[55]	76,577 *	42,001	34,576	54.85	45.15
	Average	60,011	32,089	27,923	53.12	46.88
	RMSEm	60,928.23	32,820.33	28,200.35	53.22	46.97
SVM	[56]	90,081	50,393	39,688	55.94	44.06
SAM	[57]	100,627	51,641	48,986	51.32	48.68
NNC	[46]	109,481	69,953	39,528	63.90	36.10
MXL	[58]	106,107	46,421	59,685	43.75	56.25
	Average	101,574	54,602	46,972	53.73	46.27
	RMSEm	101,839.54	55,350.20	47,695.72	54.22	46.85
	Total Average	57,484.81	30,989.19	26,495.69	54.56	45.44
	Total RMSEm	64,443.86	34,925.60	30,074.62	54.96	45.93

). Sample data used in statistical calculations for the customized NDBI approach are shown in Table 4. Average error or average bias and percent bias were calculated to produce a chart of average bias (Table 5) values. We analysed semi-automatic FE approaches for mapping of BIRs, which were categorized as [a] the customized NDBI approach (total RMSE = 2617.90 ± 58.26 m², average RMSE = 654.48 ± 58.26 m²), [b] the TD approach (total RMSE = 3951.25 ± 55.05 m², average RMSE = 987.81 ± 55.05 m²), [c] the SP approach (total RMSE = 5308.37 ± 127.83 m², average, RMSE = 1327.09 ± 127.83 m²), and [d] the PSC approach (total RMSE = 9037.70 ± 115.36 m², average RMSE = 2259.43 ± 115.36 m²). Our experiment revealed that the new customized NDBI approach outperformed the other three existing approaches for extraction of BIRs. The performance of the TD and SP approaches were comparable. The overall reliability of all the extraction methods within the respective approaches is based on various statistical parameters and is discussed below. The variation in statistical indicators of accuracy between the four approaches is summarized in Table 5 and Table 6, and Figure 5. A sample visual analysis/comparison of the blue-ice extraction results achieved by the four feature extraction approaches is given in Figure 6 and Figure 7.

(i)

Performance of the customized NDBI approach

The NDBI-1 (RMSE = 491.65 ± 58.26 m²) combination yielded superior results, while NDBI-4 (RMSE = 768.46 ± 58.26 m²) produced inferior results compared to the remaining two combinations in the four analysed NDBIs (Table 5, RMSE). Conversely, the differences in average RMSE for NDBI-2 (682.58 ± 58.26 m²) and NDBI-3 (675.21 ± 58.26 m²) were comparatively negligible in a given group. The customized NDBI approach was superior to the existing three FE approaches. It is also evident that the NDBIs (NDBI-1 and NDBI-2) designed using the green band performed slightly better than those using the yellow band (NDBI-3 and NDBI-4) (Figure 6). Similarly, NDBIs designed using the NIR-1 (NDBI-1 and NDBI-3) band performed better than those designed using the NIR-2 (NDBI-2 and NDBI-4) band. This implies that the normalization of green and yellow bands against NIR-1 yielded superior results than normalizing these spectral bands against NIR-2. Interestingly, the underestimation represented by positive bias (average 0.21 ± 0.01%, total 0.83 ± 0.01%) has contributed more in comparison to overestimation represented by negative bias (average 0.15 ± 0.03%, total 0.60 ± 0.03%) for misclassification of blue ice area (total bias) caused by all the four NDBI methods. In terms of misclassified number of pixels, the NDBI approach rendered a total of 58.02 ± 4.85% (average 59.07 ± 4.85%) overestimated versus 41.98 ± 4.85% (average 40.93 ± 4.85%) underestimated pixels out of the total misclassified number of pixels. Overall, the customized NDBI approach surpassed the accuracy of the existing three approaches.

(ii)

Performance of the target detection (TD) approach

Amongst the four TD methods applied to the WV-2 images, the MT-TCIMF (RMSE = 860.71 ± 55.05 m²) outperformed the other methods, while CEM (RMSE = 1113.45 ± 55.05 m²) performed the worst in this cohort of methods for extracting BIRs (Table 5, RMSE). ACE (RMSE = 941.54 ± 55.05 m²) delivered better results than OSP (RMSE = 1035.35 ± 55.05 m²). The underestimation represented by positive bias (average 0.32 ± 0.03%, total 1.27 ± 0.03%) contributed more in comparison to overestimation represented by negative bias (average 0.29 ± 0.02%, total 1.15 ± 0.02%) for misclassification of the blue-ice area (total bias) caused by all four TD methods.

MT-TCIMF (positive bias 0.25 ± 0.03%, negative bias 0.27 ± 0.02%) and OSP (positive bias 0.31 ± 0.03%, negative bias 0.32 ± 0.02%) caused more overestimation compared to underestimation, while ACE (positive bias 0.33 ± 0.03%, negative bias 0.24 ± 0.02%) and CEM (positive bias 0.39 ± 0.03%, negative bias 0.32 ± 0.02%) caused more underestimation compared to overestimation. In terms of misclassified number of pixels, the TD approach rendered a total of 52.46 ± 2.33% (average 52.31 ± 2.33%) overestimated versus 47.54 ± 2.33% (average 47.69 ± 2.33%) underestimated pixels out of the total misclassified number of pixels (Figure 6). In general, the performance of the TD approach in terms of average RMSE was superior in comparison to the SP and PSC approaches, and inferior to the NDBI approach.

(iii)

Performance of the spectral processing (SP) approach

The MTMF spectral processing method (RMSE = 1119.16 ± 127.83 m²) produced the best results when compared with the other methods within this approach. The MF/SAM (RMSE = 1152.36 ± 127.83 m²) ratio method performed better than its individual components, i.e., SAM (RMSE = 2179.90 ± 127.83 m²) and MF (RMSE = 1361.74 ± 127.83 m²), as expected (Table 5, RMSE). The difference in RMSE between MF/SAM and MF was comparable to the difference in RMSE between MF/SAM and SAM, which suggests that the poor performance of MF/SAM may be attributed to the overall poor performance of the SAM component. The underestimation represented by positive bias (average 0.45 ± 0.06%, total 1.81 ± 0.06%) has contributed more in comparison to overestimation represented by negative bias (average 0.39 ± 0.03%, total 1.57 ± 0.03%) for misclassification of blue ice area (total bias) caused by all four SP methods (Figure 7). It is pertinent to note that the differences in average RMSE for MTMF (1119.16 ± 127.83 m²) and MF/SAM (1152.36 ± 127.83 m²) were negligible in this group of methods. Interestingly, PCA performed significantly less well than MF/SAM, MF and MTMF. Notably, the MTMF (positive bias 0.34 ± 0.06%, negative bias 0.37 ± 0.03%) caused more overestimation compared to underestimation, while other three methods in the cohort (MF/SAM, MF and PCA) caused more underestimation in comparison to overestimation. In terms of misclassified number of pixels, the SP approach resulted in a total 53.47 ± 1.77% (average 53.12 ± 1.77%) overestimated versus 46.53 ± 1.77% (average 46.88 ± 1.77%) underestimated pixels out of the total misclassified number of pixels. In general, the performance of the SP approach in terms of average RMSE was superior in comparison to the PSC approach and inferior to the TD and NDBI approaches.

(iv)

Performance of the pixel-wise supervised classification (PSC) approach

The SVM (RMSE = 1976.81 ± 115.36 m²) combination yielded superior results, while NNC (RMSE = 2503.42 ± 115.36 m²) performed least well, compared to the remaining two combinations in this group of four pixel-based classification methods (Table 5, RMSE). Conversely, the differences in average RMSE for SAM (2179.90 ± 115.36 m²) and MXL (2377.57 ± 115.36 m²) were negligible in this group. The underestimation represented by positive bias (average 0.77 ± 0.07%, total 3.08 ± 0.07%) contributed more in comparison to overestimation represented by negative bias (average 0.66 ± 0.07%, total 2.65 ± 0.07%) for misclassification of blue ice area (total bias) of all four PSC methods (Figure 7). Notably, the MXL (positive bias 0.65 ± 0.07%, negative bias 0.84 ± 0.07%) caused more overestimation compared to underestimation, while the other three methods in the cohort (SVM, SAM, NNC) caused more underestimation in comparison to overestimation. In terms of misclassified number of pixels, the PSC approach resulted in a total of 53.76 ± 4.22% (average 53.73 ± 4.22%) overestimated versus 46.24 ± 4.22% (average 46.27 ± 4.22%) underestimated pixels out of the total misclassified number of pixels. The PSC approach was inferior to the other three FE approaches (SP, TD, NDBI) analysed here.

(v)

Overall performance of semi-automatic extraction methods

The overall performance for all the blue ice extraction methods, based on RMSE and bias, is summarized in Table 5, Table 6, and Figure 5 and was ranked as follows: NDBI-1 > NDBI-3 > NDBI-2 > NDBI-4 > MT-TCIMF > ACE > OSP > CEM > MTMF > MF/SAM > MF > PCA > SVM > SAM > MXL > NNC. This order suggests that the methods can be grouped by approach and ranked as follows: NDBI approach > TD approach > SP approach > PSC approach. The RMSE values for the customized NDBI approach ranged from ~491 ± 58 to ~769 ± 58 m², and for the TD, SP and PSC approaches, RMSE values varied from ~860 ± 55 to 1114 ± 55 m², ~1119 ± 128 to ~1675 ± 128 m², and ~1976 ± 115 to ~2503 ± 115 m², respectively (Table 5, Figure 5). The percentage bias values for the customized NDBI approach ranged from ~0.28 ± 0.03 to ~0.40 ± 0.03%, and for the TD, SP and PSC approaches, from ~0.52 ± 0.04 to ~0.71 ± 0.04%, ~0.71 ± 0.09 to ~1.08 ± 0.09%, and ~1.27 ± 0.06 to ~1.54 ± 0.06%, respectively. The smallest variation in RMSE and bias was found in the customized NDBI approach, which suggests that it is a more stable and consistent approach than the other three blue ice extraction approaches. In terms of the internal stability associated with each approach, they ranked as follows: customized NDBI > TD > SP > PSC. Along with the NDBI methods, the MTMF methods (MTMF of SP approach and MT-TCIMF of TD approach) outperformed the other extraction methods in their respective cohorts, while the CEM, PCA and NNC methods had the poorest performance of the 16 extraction methods. In general, NDBI methods consistently achieved a minimum number of false positives for all 12 tiles.

Our results support the strong potential of the NDBI approach for extracting BIRs. Among all the TD, SP, and PSC methods, excluding MTMF methods, results of the ACE method showed great target visibility and the capability to suppress false signals, resulting in a low RMSE (941.54 ± 55.05 m²) when compared with other methods. The comparative analysis of RMSEs (Table 5) suggests that MT-TCIMF (TD approach) (RMSE = 860.71 ± 55.05 m²) was able to greatly reduce the number of false signals when compared with MTMF (SP approach) (RMSE = 1119.16 ± 127.83 m²). Our results are consistent with earlier studies, which revealed that if the spectral angle between the target and the non-target is significant, MT-TCIMF may reduce the number of false positives when compared with the CEM method. Based on the RMSEs computed from the number of false assignments, ACE and MT-TCIMF achieved better performance than the SP and PSC approaches. We also found that the TD and SP methods were suitable for blue ice detection, even though their performance was poorer than that of NDBI methods. The MF/SAM ratio suppressed the false positives that were present when using only one of the methods (MF or SAM) but not when using the other, while enhancing the true positives. For example, if a pixel representing blue ice had a high MF and low SAM value, the ratio (MF/SAM) would produce a high true positive. Conversely, if the MF product had a high value for a false positive, but the SAM correctly mapped it as a non-blue ice pixel (high value), the ratio of these high values led to a smaller value, which suppressed the false positive. Hence, MF/SAM (RMSE = 1152.36 ± 127.83 m²) performed considerably better than MF (RMSE = 1361.74 ± 127.83 m²) and SAM (RMSE = 2179.90 ± 115.36 m²) individually.

4.1.2. Topographical Influences on the Extraction of BIRs

It is assumed that the use of normalised band ratios considerably reduces the topographical effect in comparison to the use of single spectral bands [59]. However, a few studies have shown that the topography significantly affects band ratios in hilly areas [60,61,62] and image classification [63]. Therefore, it is pertinent to demonstrate the robustness and validity of the analysis by testing the effect of topography. We quantitatively analysed the effects of the topography of the study area on the extraction of BIRs using four feature extraction approaches. The topographic distribution of the 12 tiles analysed is highlighted in Table 2. To test the robustness of the analysis, we describe the effect of topography based on three parameters.

(i)

Comparative analysis of errors for various tiles in different elevation settings

RMSE values for three types (LET, MET, and HET) of tiles (based on elevation) and local statistics for four semi-automatic FE approaches are given in Supplementary Table S3. Variation of the percent bias, including overestimated area (negative bias) and underestimated area (positive bias), was estimated. The contribution of overestimated pixels and underestimated pixels (%) to the total bias in the extracted blue-ice area is shown in Figure 8. To test the performance of the four FE approaches in blue ice extraction in terms of overall stability and consistency, box plots were constructed (Supplementary Figure S6).

The average and total RMSE values for each FE approach (column) were calculated to assess their overall stability and consistency (Supplementary Table S3). The 16 individual FE methods were categorized into four approaches as (a) the customized NDBI approach (average RMSE, MET = 562.72 ± 95.45 m², HET = 264.78 ± 51.18 m², LET = 1015.23 ± 85.83 m²), (b) the TD approach (average RMSE, MET = 868.66 ± 73.38 m², HET = 666.49 ± 128.57 m², LET = 1441.51 ± 57.84 m²), (c) the SP approach (average RMSE, MET = 1182.64 ± 67.32 m², HET = 1048.44 ± 229.17 m², LET = 1807.61 ± 124.62 m²), and (d) the PSC approach (average RMSE, MET = 1978.97 ± 187.32 m², HET = 1881.32 ± 224.80 m², LET = 2978.58 ± 222.60 m²). Our analyses revealed that the new customized NDBI approach outperformed the three existing approaches for extraction of blue ice areas from various types of elevated areas. This implies that the NDBI approach is least prone to topographic and illumination-based errors in classifying blue ice areas. The TD and SP approaches showed comparable though less effective performance. It was also evident that the NDBI approach led to the least errors for HET (average RMSE = 264.78 ± 51.18) as compared to MET (average RMSE = 562.72 ± 95.45) and LET (average RMSE = 1015.23 ± 85.83), indicating that the NDBI approach is relatively less precise in extracting blue ice from LET compared to MET and HET.

The overall trends of performance for all the blue ice extraction methods, based on RMSE and bias, are summarized in Figure 8 and Figure S6. It is notable that the accuracy trend for extracting blue ice from topographically varying tiles is consistent for all four approaches, in the order HET > MET > LET. This implies that the 16 pixel-based methods were more suitable to extract BIRs from HET compared to LET and MET, which can be attributed to the sun position and surface illumination at the time of WV-2 image acquisition. It is also evident that the different topographical tiles rendered varying blue ice mapping accuracies for the different methods. The average RMSE values for the customized NDBI approach for all three topographical categories (LET, MET, HET) ranged from ~264 ± 77 to ~1015 ± 77 m², and for the TD, SP and PSC approaches, average RMSE values varied from 666 ± 87 to 1441 ± 87 m², from ~1048 ± 140 to ~1807 ± 140 m², and from ~1881 ± 22 to ~2978 ± 21 m², respectively (Supplementary Table S3). The smallest variation in RMSE and bias was found in the customized NDBI approach, which suggests that it is the most stable and consistent approach. In terms of the internal stability associated with each approach, we rank these as follows: customized NDBI > TD > SP > PSC. Figure 8 shows the variation of percent bias values including positive bias and negative bias for extracting blue ice from topographically varying tiles (LET, HET, MET) using the 16 pixel-based methods categorized in four FE approaches. Having concluded above that HET yielded comparatively higher accuracy for all the 16 methods, it is interesting to note that misclassification in all the 16 methods was almost completely because of overestimation (negative bias) rather than underestimation (accounting for almost 0%). Misclassification in MET was mostly caused by underestimated pixels (positive bias), while that in LET was mostly caused by overestimated pixels (negative bias).

(ii)

Comparative analysis of errors for shadow-covered against non-shadow covered tiles.

Of the 12 tiles under consideration, two (tiles 8, 9) were in shadow-prone regions as identified by digital elevation model (DEM)-based hill shade (shaded relief) analysis. Supplementary Figure S7 presents the RMSE results pertaining to those two shadowed tiles (ST) and 10 non-shadowed tiles (NST) (tiles 1, 2, 3, 4, 5, 6, 7, 10, 11, 12). The overall trend of accuracy for the blue ice extraction methods as they pertain to shadow cover, based on RMSE, was: NDBI-1 > NDBI-3 > NDBI-4 > NDBI-2 > MT-TCIMF > ACE > OSP > CEM > MTMF > MF/SAM > MF > PCA > SVM > SAM > MXL > NNC. This ranking suggests that the methods can be grouped by approach and ranked based on their ability to extract shadow-covered blue ice is as follows: NDBI > TD > SP > PSC (Supplementary Figure S7). For shadowed tiles, the average RMSE for the customized NDBI approach varied from ~946 ± 117 to ~1382 ± 117 m². The RMSE of the TD, SP and PSC approaches varied from ~1418 ± 61 to ~1715 ± 61 m², ~1691 ± 155 to ~2394 ± 155 m², and 2866 ± 323 to ~4250 ± 323 m², respectively. The overall trend of accuracy for the blue ice extraction methods for non-shadowed tiles, based on RMSE, was: NDBI-1 > NDBI-2 > NDBI-3 > NDBI-4 > MT-TCIMF > ACE > OSP > CEM > MTMF > MF/SAM > MF > PCA > SVM > NNC > SAM > MXL. This ranking suggests that the methods can be grouped by approach based on their ability to extract blue ice from NST as follows: NDBI > TD > SP > PSC approach (Supplementary Figure S7). For non-shadowed tiles, the average RMSE for the customized NDBI approach varied from ~333 ± 100 to ~784 ± 100 m², which is significantly lower than that for the shadowed tiles (~946 ± 117 to ~1382 ± 117 m²). The RMSE of the TD, SP and PSC approaches varied from ~697 ± 58 to ~948 ± 58 m², ~964 ± 124 to ~1489 ± 124 m², and ~1745 ± 83 to ~2144 ± 83 m², respectively. These observations clearly confirm that the shadowed regions cause more misclassification than the non-shadowed regions in the extraction of blue ice areas.

(iii)

Comparative analysis of errors based on amount of blue ice area present in various tiles

The spatial distribution of tiles based on surface area of BIRs is described in Table 2, categorized into three groups: LBT (comparatively low blue ice, tiles 6, 8, 9), HBT (relatively high blue ice, tiles 1, 3, 4, 5, 7), MBT (intermediate blue ice, tiles 2, 10, 11, 12). Statistics were calculated separately for these three categories to compare the influence of amount of blue ice in different tiles (Supplementary Table S4). Supplementary Figure S8 shows the variation of RMSE rendered by 16 FE methods categorized into four approaches for extraction of blue ice from various tiles based on their classification as LBT, MBT or HBT). The average and total RMSE values for each FE approach (column) was calculated to assess their overall stability and consistency (Supplementary Table S4). Our analyses revealed that the new customized NDBI approach outperformed the three existing approaches for extraction of blue-ice areas from tiles varying in the amount of blue ice included. Supplementary Figure S8 shows that tiles with moderate or high blue-ice percentages give comparable accuracy for all 16 extraction methods, while tiles with low blue ice percentage are prone to misclassification, achieving lower accuracy in terms of RMSE. These observations clearly confirm that the accuracy of pixel-based methods is critically dependent on the amount of the target feature and confounding non-target features present in the study area.

4.2. Comparing Areas of Extracted BIRs against Reference BIRs for the Entire Study Region

In the previous section, we demonstrated the robustness of the analysis using the data collected on the field using differential GPS supported by manual digitization. In order to compare the blue ice extraction results with existing blue ice maps and to test the validity against solely manually digitised reference map, we compared our extraction results against the entire extent of the study. A pre-digitized database of the blue ice areas in the entire study region, in the form of shapefiles, was utilized as a reference for assessing the accuracy of the semi-automatic extraction methods. Using the vectorized outputs from the extraction procedures, the areas of the extracted blue ice areas were estimated, and area statistics were compared with blue ice areas in the manually digitized database. The reference data used for this analysis are given in

Table 7. Manually digitized data summarizing the reference blue ice regions (BIRs) and non-BIRs (snow/supraglacial debris, streams, etc.) for entire study area. The reference blue-ice area derived from existing MODIS/Landsat ETM+ [22] data is also described.

Reference Data		Reference Blue Ice Data
Reference Map	Total Area (m2)	BIRs (m2)	Non-BIRs (m2)	% BIRs	% Non-BIRs
Manually digitized BIRs	173,530,259.00	106,875,250.60	66,655,008.40	61.59	38.41
MODIS–ETM+ BIRs	173,530,259.00	140,281,789.02	33,248,469.98	80.84	19.16

.

The extracted BIR map was compared with manually-digitized reference data to produce a chart of total bias (m² and %) (

Table 8. Quantitative evaluation (bias, % and m²) of various feature extraction methods for extracting blue ice from the entire study area. The lowest (italics) and highest (bold) values in each local column are highlighted. The column-wise average values are bolded and underlined. Total bias represents the total misclassified blue ice area in % and m². Local outliers are highlighted using *.

		Extracted Blue Ice Area	Total Bias (m2)	Total Bias (%)	Remark
NDBI-1	Present work	106,253,876.47 *	621,374.13 *	0.58 *	Underestimation
NDBI-2	Present work	107,510,771.62	−635,521.02	−0.59	Overestimation
NDBI-3	Present work	107,493,642.98	−618,392.38	−0.58	Overestimation
NDBI-4	Present work	107,572,503.29	−697,252.69	−0.65	Overestimation
	Average	107,207,698.59	643,135.06	0.60
	RMSEm	107,209,116.94	643,926.09	0.60
MT-TCIMF	[48]	108,050,194.52	−1,174,943.92	−1.10	Overestimation
OSP	[49]	108,069,168.73	−1,193,918.13	−1.12	Overestimation
ACE	[50]	105,606,034.11	1,269,216.49	1.19	Underestimation
CEM	[51]	105,584,361.53	1,290,889.07	1.21	Underestimation
	Average	106,827,439.72	1,232,241.90	1.15
	RMSEm	106,834,546.85	1,233,211.13	1.16
MTMF	[52]	108,595,298.82	−1,720,048.22	−1.61	Overestimation
MF/SAM	[53]	108,639,032.38	−1,763,781.78	−1.65	Overestimation
MF	[54]	105,082,210.18	1,793,040.42	1.68	Underestimation
PCA	[55]	105,013,550.34	1,861,700.26	1.74	Underestimation
	Average	106,832,522.93	1,784,642.67	1.67
	RMSEm	106,847,432.04	1,785,386.08	1.67
SVM	[56]	104,108,595.61	2,766,654.99	2.59	Underestimation
SAM	[57]	104,041,038.89	2,834,211.71	2.65	Underestimation
NNC	[46]	103,906,014.77	2,969,235.83	2.78	Underestimation
MXL	[58]	109,840,008.48 *	−2,964,757.88 *	−2.77 *	Overestimation
	Average	105,473,914.44	2,883,715.10	2.70
	RMSEm	105,504,057.77	2,885,016.75	2.70
	Total Average	106,585,393.92	1,635,933.68	1.53
	Total RMSEm	106,600,768.17	1,833,465.00	1.72

). Considering the customized NDBI approach, NDBI-1 (total bias = 0.58 ± 0.02%) and NDBI-3 (total bias = −0.58 ± 0.02%) yielded superior results, while NDBI-2 (total bias = −0.59 ± 0.02%) and NDBI-4 (total bias = −0.65 ± 0.02%) produced inferior results in the group of four analysed NDBIs (Table 8, Total bias). The differences in total bias for NDBI-1 (0.58 ± 0.02%), NDBI-2 (−0.59 ± 0.02%) and NDBI-3 (−0.58 ± 0.02) were comparatively negligible in a this group. However, NDBI-1 resulted in underestimation of blue ice, while the remaining three NDBI models resulted in overestimation. The NDBIs (NDBI-1 and NDBI-2) designed using the green band performed slightly better than those designed using the yellow band (NDBI-3 and NDBI-4). Similarly, NDBIs designed using the NIR-1 band (NDBI-1 and NDBI-3) performed better than those designed using the NIR-2 band (NDBI-2 and NDBI-4). Overall, the customized NDBI approach surpassed the accuracy associated with the existing three approaches, suggesting that customized NDBI resulted in fewer false assignments. Considering the TD approach, MT-TCIMF (total bias = −1.10 ± 0.03%) and OSP (total bias = −1.12 ± 0.03%) yielded superior results, while ACE (total bias = 1.19 ± 0.03%) and CEM (total bias = 1.21 ± 0.03%) produced inferior results, in this group of four analysed methods (Table 8, Total bias). The differences in total bias for MT-TCIMF/OSP and ACE/CEM pairs were also comparatively negligible in this group. It is notable that MT-TCIMF and OSP resulted in overestimation of blue ice, while ACE and CEM resulted in underestimation. Considering the SP approach, MTMF (total bias = −1.61 ± 0.03%) and MF/SAM (total bias = −1.65 ± 0.03%) yielded superior results, while MF (total bias = 1.68 ± 0.03%) and PCA (total bias = 1.74 ± 0.03%) produced inferior results, in this group of four analysed methods (Table 8, Total bias). The differences in total bias for MTMF/MF/SAM pair were again comparatively negligible in this group. In this group, MTMF and MF/SAM resulted in overestimation of blue ice, while MF and PCA resulted in underestimation. Finally, considering the PSC approach, SVM (total bias = 2.59 ± 0.05%) and SAM (total bias = 2.65 ± 0.05%) yielded superior results, while MXL (total bias = −2.77 ± 0.05%) and NNC (total bias = 2.78 ± 0.05%) produced inferior results, in this group of four analysed methods (Table 8, Total bias). The differences in total bias for the NNC/MXL pair were also comparatively negligible in this group. In this group, MXL resulted in overestimation of blue ice, while the other three methods in this cohort resulted in underestimation.

The overall trend of accuracy for the blue ice extraction methods, based on total bias (%), was: (NDBI-1 = NDBI-3) > NDBI-2 > NDBI-4 > MT-TCIMF > OSP > ACE > CEM > MTMF > MF/SAM > MF > PCA > SVM > SAM > MXL > NNC. This analysis is very consistent with the tile-wise analysis described in the previous section. A visual image analysis of extracted BIRs from the entire study area using the NDBI approach is shown in Figure 9 (TD approach on Supplementary Figure S9, SP approach on Supplementary Figure S10 and PSC approach on Supplementary Figure S11). The visual analysis reveals that most instances of misclassification occurred at low elevation or in shadow prone areas. Most of the feature extraction methods failed to recognize shadowed blue ice pixels. However, it was notable that the NDBI-derived shadowed blue ice pixels were relatively less misclassified, exemplifying the robustness of the NDBI approach for extraction of blue ice in shadowed areas. In addition, most misclassification was constrained to areas of dense sastrugi, which are sharp irregular grooves/ridges formed on a snow surface by wind erosion, saltation of snow particles and deposition. These regular ridges cause local shadowing because of rough topography. These local patches of sastrugi hamper the detection of blue ice pixels. The blue ice area related to seasonal melting and ice flows (“rough BIRs”) has a relatively rough surface, with troughs and ridges caused by melting and water flows, and crevassing caused by high strain in the ice. The troughs and crevasses trap snow and lead to mixed pixels that contain both blue ice and snow. The ridges also induce a solar shadow effect on the WV-2 imagery. In addition, wind crusts present micro-relief with relatively smaller patches mixed with snow, again forming a rough surface texture, which was the major source of misclassification in the present study.

4.3. Comparison of the Extracted BIRs with Existing BIR Map

In the final step, semi-automatically extracted blue ice areas from WV-2 were compared with the existing blue ice map derived from MODIS and ETM+ satellite data [22] (Figure 10). It is evident that the MODIS and ETM+ derived BIR map overestimated the blue ice in the study area. Considering the coarse resolution of MODIS and medium resolution of ETM+, the undersized snow cover patches were misclassified as blue ice. This necessitated refinement of the existing BIR map derived from medium resolution satellite data by VHR satellite data. The present study aimed to refine the existing local BIR map using VHR satellite data to exclude the local snow cover. It is also evident that the MODIS and ETM+ derived BIR map includes some of the lakes of Schirmacher oasis as blue ice, which needs to be rectified. Epishelf lakes and frozen lakes, in particular, confound the spectral signature of blue ice leading to their misclassification as blue ice. The WV-2 derived map clearly excludes the local snow cover and lakes of the Schirmacher Oasis, and hence leads to improvement in the local BIR map. In Figure 10, regions of overestimation of BIRs (by snow cover) derived from MODIS and ETM+ are clearly visible. These regions of overestimation are refined in the BIR map derived from WV-2 data and quantitatively, the WV-2 derived BIR map refined the existing MODIS/ETM+ derived BIR map by around 31% (~33,406,538 m²).

4.4. Errors Associated with Manually Digitized Reference Data

Assessment of the accuracy of the digitized vector data was necessary as this is the most crucial reference data used in this study. A perpendicular error distance from the digitized shoreline vertex to the DGPS reference shoreline vertex was calculated for 10 random points per blue ice tile, i.e., 120 points from 12 tiles. The perpendicular distance errors varied from ~3 mm to ~11 mm (RMSE ~9 mm). Additionally, to check the effect of digitization errors on overall blue ice area, each digitized blue ice polygon area value was compared with the original reference DGPS-based blue ice area value. The difference in blue ice area values varied from ~4 m² to ~7 m² with RMSE of ~8.5 m². These quantitative errors in digitization are insignificant in comparison to the errors resulting from the different extraction methods (~200 m² to ~3000 m²) for the extraction of BIRs. Furthermore, the digitization error would be constant across all 16 methods and, hence, comparison of these methods would be unbiased and valid.

5. Discussion

5.1. Testing the Performance of Supervised Information Extraction Methods in the Antarctic Environment

Overall, the visual images of the extracted BIRs obtained using four general feature extraction approaches indicated that most misclassification occurred in shadow-prone and topographically varying areas (Figure 8 and Figures S6–S11). We also noted positive (representing exceptional accuracy in BIR mapping in this cohort) and negative (representing inferior performance in this cohort) outliers. These outliers are highlighted in Table 4, Table 5, Table 6, Table 8, Tables S3 and S4, and related to all four approaches. Global outliers (w.r.t. the 16 individual methods) and local outliers (w.r.t. the four methods in each approach) are highlighted in Table 6. The occurrence of outliers was related to excessive overestimation or excessive underestimation of blue ice area. Outliers found in the present study are significant as they represent the factors affecting the performance of different methods. For instance, NDBI-4 yielded a significant underestimation (compared to all 12 tiles) for deriving BIRs from shadow-prone tiles (tile 8 and tile 9), indicating its poor ability to map BIRs with many shadows. NDBI-1 provided much superior performance in a cohort of four models representing the positive outlier. The exceptional overestimations in MXL appear to be caused by confounding signatures of background features. Poor performance and misclassifications by MXL have been reported previously [64,65,66]. This supports the view that non-parametric methods such as SVM are better suited to VHR compared to parametric methods such as MXL [67]. The NDBI approach used in the present study did not suffer from significant negative outliers and hence performed consistently throughout the experiment. Our analyses also confirm previous studies [45] which have suggested that MT-TCIMF can significantly reduce the number of false positives when compared with CEM.

5.2. Effect of New Spectral Bands of the WV-2 Used in NDBI for Blue Ice Mapping

The NIR-2 and yellow band of WV-2 are the strategic bands used to design NDBIs, and are not available through any other satellite platforms. This provides a significant advance through the new spectral bands provided by WV-2. The four NDBIs trialled here can potentially be differentiated into two traditional NDBIs using the conventional Green and NIR-1 bands and two new NDBIs generated using the yellow and NIR-2 bands. Our analyses show that the green and NIR-1 bands provide comparable results to the new spectral bands. Therefore, contrary to several studies that have reported significant improvements in mapping various targets using the new spectral bands of WV-2 (e.g., [68,69,70,71,72]), our study did not identify any significant improvement in mapping blue ice using these new spectral bands.

5.3. Comparison with Previous Blue Ice Mapping in Antarctica

To date, there have been 16 dedicated studies [2,15,16,21,22,23,25,29,30,31,33,34,35,42,73,74] published between 1990 and 2022 (seven before and nine after the launch of WV-2) which focus on blue ice mapping applications. Of these, only five [2,22,24,35,73] attempted a continent-wide approach. All other studies mapped BIRs in limited local areas. None of these studies used VHR data. Of the 16 studies, only one [34] was thoroughly validated using ground-truth data while the remainder included qualitative analysis of results based on limited ground survey observations and visual interpretations. The present study’s results were validated using qualitative (visual interpretation) and thorough quantitative (ground-truthing) analyses. Most of the previous studies were based on use of one or a combination of two methods, including optical grain-size thresholding, SIRs, unsupervised classification, object-oriented classification, thresholding albedo values and MXL supervised classification. The current study is the first systematic study to test the performance of 16 BIR mapping methods quantitatively and qualitatively in a significant region in Antarctica. Previous case studies using MODIS and Landsat yielded errors ranging from 2–4 km² [15], while an ASTER-based study yielded errors ranging from 0.3–0.7 km² [34]. In contrast, our study yielded errors of 0.0049 km² to 0.0274 km². This is consistent with the study of Yu et. al. [30], who reported the improved lower estimates of BIR extent in the Lambert Glacier–Amery Iceshelf system using Landsat ETM+ as compared to those derived by using AVHRR.

During the summer season, daily BIR maps can be generated using a coarse-resolution MODIS images, while weekly or monthly maps can be generated using Landsat and Sentinel. However, precise and refined mapping of BIRs can only be achieved by using high spatial-spectral resolution satellite imageries, as demonstrated in the present study. An alternative strategy could involve mapping larger BIR patches using MODIS, moderate-sized patches using Landsat/Sentinel and small-sized patches using VHR data. For a continent-wide mapping, Landsat and MODIS provide a good choice. However, to study regional and sub-regional dynamics with high accuracy, WV-2 based mapping is appropriate. The maximum summer extent of BIRs should be estimated using WV-2 across Antarctica.

5.4. Statistical Significance and Generalization of Performance of Methods

The extent of misclassification arising from the 16 methods in terms of bias in the extracted blue-ice areas varied from 4,931 m² to 27,370 m², with differences in bias ranging from 0.1% to 1.26%. Noting that a 1% bias (~1652 m²) in area caused a misclassification of ~17,734.2 m² (~70,937 pixels) blue ice, even small changes in bias values represent considerable errors in misclassification. The use of high spatial resolution WV-2 data inherently leads to an expectation of improved accuracy in the mapping of blue ice. In our analyses, the small changes in bias values (0.1 to 0.9%) are therefore significant in the mapping of BIRs. The mathematical ranking of methods based on small variations in bias values is statistically stable and practically robust. Considering the most recent estimate of the total area of BIRs in Antarctica of 234,549 km² (derived by Hui et al. [22] using LIMA), if total BIR were to be derived across Antarctica using WV-2 and upscaling up our methodology, the errors in blue-ice areas would vary between 234.5 km² and 3612.0 km² using the 16 methods.

5.5. Beyond Quantitative Analysis and Evaluation of Accuracies

A methodology that is robust, accurate and user-friendly can reduce time-consuming and laborious manual digitization. However, a semi-automatic blue-ice FE can be applied in a real operational environment only if it provides superior performance in terms of quality measures. In addition to empirical superiority in terms of accuracy, our analyses showed that the customized NDBI provides the following advantages over existing extraction methods.

(1)

Consistency: customized NDBI-based blue-ice mapping was less sensitive to the background noise (topography, shadows, snow/ice cover) and could consistently extract blue ice from the study area by optimizing the target against the noise using green, yellow, NIR-1 and NIR-2 bands;

(2)

Flexibility: the slight variation of the input approximations (threshold values in the case of customized NDBI and ROIs in cases of SP, PSC or TD approaches) should not affect the extraction results significantly. The threshold values used for customized NDBI can be adjusted to suppress the noise with no effect on the extraction of targeted blue ice;

(3)

Minimizing of manual editing: a semi-automatic method should eventually minimize manual editing of the extracted BIRs. A visual interpretation of the extracted blue ice showed that the dimensions of BIRs (shape, geometry and size) were well-preserved and manual editing was effectively minimized;

(4)

Efficiency: the extraction should be executed, and the results should be available rapidly. The total time required for the extraction should be much less than that by manual digitization. Efficiency is also dependent on reliability, accuracy, and interactivity factors. The average time for extracting the blue ice was the lowest for the customized NDBI, moderate for SP and PSC, and maximal for TD approaches;

(5)

Interactivity: the semi-automatic extraction strategy should be an interactive process between the machine learning algorithms and the human operator. This should allow the operator to correct erroneous results immediately after extraction. The customized NDBI approach is highly interactive in terms of threshold definition, and in band selection to rectify wrongly classified pixels;

(6)

Robustness: the semi-automatic blue ice mapping method should work well for different types of BIRs in the cryospheric environment. Most of the BIRs in the study area here were influenced by local topography. The customized NDBI method could effectively identify BIRs in various topographical settings in this typical Antarctic coastal oasis environment;

(7)

Complexity: any semi-automatic method should ideally be simple to implement on satellite images. The customized NDBI approach is simpler to implement than either the TD or SP algorithms;

(8)

Accuracy: extracted information should be correct and geometric errors should be minimized, giving results that are at least comparable to those achieved by manual digitization. The RMSE values obtained from the customized NDBI confirmed that the approach could extract blue ice accurately w.r.t. the manually digitized reference;

(9)

Visual comparison: the semi-automatic method should provide extraction results which can be compared visually against the manual reference data. Visual comparison here demonstrated that blue-ice polygons from all 12 tiles were detected in the NDBI images, and that the boundaries of the extracted blue ice polygons matched the actual boundaries of the blue ice in the images or reference digitized data closely;

(10)

Error: variation in error should be minimized. A box plot (Supplementary Figures S6 and S7) shows that the customized NDBI had superior performance over the other three approaches tested for extracting BIRs;

In summary, our objective of developing a consistently reliable customized NDBI approach to automate previous needs for manual digitization for blue-ice mapping was achieved.

5.6. Factors Affecting the Performance of NDBI

Owing to low reflectance in the blue band and extremely low reflectance of longer wavelengths, ice present on lakes or ponds appears deep blue, with an at-satellite reflectance of 40–60% for the blue band and 20–40% for the NIR band. In contrast, BIRs show as light blue, with an at-satellite reflectance of 90–95% for the blue band and 30–40% for NIR band. The reflectance values of frozen ice on lakes or ponds are significantly lower than BIRs for all spectral bands, relatively in the NIR band, enabling the customized NDBI to differentiate blue ice from ice cover on lakes. Additionally, errors arising from shadowing were minimal in the final stage of extraction using the customized NDBI. Some manual correction was still required at the final stage, specifically including some small, shadowed patches of BIRs in mountainous and “rough BIRs” in crevassed/sastrugi/ripple areas.

The range of threshold values for NDBIs was chosen to match visual interpretations and spectral signatures of blue ice at different locations considering various confounding factors. Because of robust choice of threshold values, the NDBI was relatively insensitive to small shifts of the chosen threshold values. Overall, the NDBI method was more conservative in discriminating darker blue ice from shaded or dimly lit snow or rock–snow mixtures. However, it was somewhat sensitive to mixtures of snow over blue ice: patchy snow cover at the edges of blue ice patches (typically more than 40%) yielded reflectance values overlapping with the allowed range for blue ice. In this circumstance, our analysis showed that the NDBI provided better results than the other methods trialled.

5.7. Factors Affecting Blue Ice Mapping in Cryospheric Environments

There were seven major surface land cover or topographical factors (Figure 11) in the study area that contributed to misclassification in blue ice mapping (Supplementary Figure S12). Water bodies, supraglacial debris and streams, crevasses, and ripples were the major confounding factors causing misclassification using customized NDBI methods, while water bodies, meltwater on ice surface, and topography caused significant misclassification using TD methods (Figure 11). Performance of the SP and PSC approaches was strongly affected by all seven confounding features, making these methods more vulnerable to misclassification. The NDBI and TD approaches were not affected by landmass/rocks/nunataks, while the PSC approach was moderately impacted by this factor.

Seasonal changes in blue ice are challenging to monitor due to events such as enhanced melting, formation of melt ponds, snow events, blizzards, and even cloud cover in satellite images. The present study highlights the use of NDBI in a cloud-free satellite image obtained on one specific date. In order to achieve seasonal monitoring of blue ice using WorldView-2, repeated acquisition of satellite data will be required. Furthermore, in situ data such as meteorological data on precipitation, blizzard or storm occurrence can be used to interpret NDBI results and adjust thresholds accordingly. Such adjustments can be done case by case basis. Images before and after an event such as storm/blizzard/excess melting can be used to fill data gaps and understand the effects of such events on NDBI. Such data can be used to establish robust threshold values to be used for seasonal change studies.

5.8. Accuracy Analysis and the Measures of Uncertainties

Logistic and practical limitations mean that exhaustive field surveys to assess the accuracy of our remotely sensed BIR mapping are not possible. While manual delineation of blue ice is possible, it is very labour-intensive and time-consuming. Therefore, we focused on a 3-step accuracy analysis: (a) using 12 manually-digitized and ground-surveyed tiles; (b) using manually-digitized BIRs from the entire study area; and (c) using results from recently published literature as validation data. The BIRs derived from our study amounted to 106,875,250.60 m², which is around 20% less than the previous MODIS/ETM+-based estimate (140,281,789.02 m²) for this study area. This difference may be related to temporal changes in BIR extent, the use of different identification algorithms or the different spatial and spectral resolutions of the sensors used. The areas and spatial distributions of BIRs determined from satellite images at the same place but at different times could vary widely. For instance, smaller areas of blue ice may be exposed if strong winds result in thin snow cover while, conversely, melting of thin snow during periods of higher air temperatures could expose larger BIRs. Such changes are also likely to vary seasonally and inter-annually, although data are not available. The different spatial and spectral resolutions of the WV-2, MODIS and ETM+ data are also likely to lead to differences in BIR estimations. Pixels containing both blue ice and snow in the 250 m resolution MODIS data and the 15 m resolution Landsat ETM+ data (as well as in the 125 m gridding-space MOA data) can either be omitted or counted as BIRs, leading to under- or overestimation of BIRs. However, these pixels are clearly recognized in the 0.5 m resolution WV-2 data. Previous studies suggest that uncertainties such as mixed pixels, crevasses and shadows could account for the difference between identified BIRs [2]. In the future, 16-bit images with an increased number of spectral bands will play an important role in mapping BIRs.

5.9. Future Directions

The refined blue ice map can be used for linking science questions such as mass balance and hydrological changes in the study area. However, for this to be achieved, improved spatio-temporal datasets are required. Conducting a thorough temporal analysis of spatial extent of BIRs by compiling time series of blue ice map datasets from the existing literature was beyond the scope of this study. However, a simple comparison of our estimates of BIR extents with previously published results may highlight the effects of spatial resolution, type of sensor, and the seasonality of image acquisition. The extent of BIRs reported by Hui et al. [22] based on MODIS/Landsat ETM+ images at 250 m/15 m spatial resolution was noticeably greater than our estimate. The coarse spatial resolution of MODIS images means that some small BIRs may have been missed and shapes of BIRs could have been distorted. Although comparison with Hui et al. [22] may indicate temporal change in BIR extent, comparison of temporal GE images suggested similar sizes, and further accurate temporal studies are required to confirm the occurrence of physical changes. The above discussion suggests that spatial resolution, type of sensor, image acquisition season and classification algorithm may all contribute to differences obtained in blue ice estimates. Future reliable analyses of long-term temporal changes in BIR extent should ideally be based on image data acquired in the same season by similar sensors. At present, we conclude that WV-2 imagery provides the best available spectral information to accurately detect and map blue ice. Using the VHR WV-2 imagery, the present study was successful in precisely mapping fuzzy boundaries of BIRs. Given that our approach was based on physical principles of spectral remote sensing, it clearly has potential to be modified to other regions of the cryosphere such as Greenland, Svalbard, Alaska, etc. Furthermore, the method could be modified to map supraglacial lakes.

5.10. Transferability of Methods

NDBIs described in this study can be made fully transferable and can reproduce results by using other sensors with similar spectral band combinations. There are many sensors (Supplementary Table S6) with the availability of traditional Green and NIR-1 spectral bands similar to WV-2. However, NIR-2 and Yellow bands are not currently available in any other satellite sensor. Therefore, NDBIs generated using Green and NIR-1 are transferable to other sensors while NDBIs generated using Yellow and NIR-2 will require new satellite sensors to allow transferability.

6. Conclusions

The BIR dataset presented in this study is the first precise and refined map covering the Schirmacher Oasis and its environs using WV-2 data. The complete high-resolution blue ice distribution map derived from this study provides the baseline information required to underpin future temporal change analysis of the BIRs in this region. The total area of blue ice in the study area is estimated to be 106.875 km², ~61% of the study area extent. The WV-2 derived BIR map clearly excludes the local snow cover (windswept firn) and various types of lakes present in the Schirmacher Oasis and gives a considerable improvement in the local BIR map. Regions of overestimation of BIRs (by snow cover) derived from MODIS and ETM+ are clearly visible. These regions of overestimation are refined in the present BIR map derived from WV-2 data. Quantitatively, the WV-2 derived BIR map refined the existing MODIS/ETM+ derived BIR map by ~31% (~33,406,538 m²). The use of all customized NDBI combinations derived by using the duplet set of VNIR bands offers a rapid and precise means for extracting BIRs as compared to existing pixel-based supervised methods. Near-nadir, mid-summer viewing conditions for the WV-2 sensor should be used in future for time series WV-2 image analysis of blue ice extent using a normalized difference algorithm. The outcome of this study is an NDBI algorithm that can be used for high precision long-term monitoring of blue-ice extent changes. The detailed and accurate map of blue ice generated in this study will provide useful input to future studies of mass balance, solar radiation budget and regional climate changes in the study region.