Comparison of Three Methods for Distinguishing Glacier Zones Using Satellite SAR Data

Abstract

Changes in glacier zones (e.g., firn, superimposed ice, ice) are good indicators of glacier response to climate change. There are few studies of glacier zone detection by SAR that are focused on more than one ice body and validated by terrestrial data. This study is unique in terms of the dataset collected—four C- and L-band quad-pol satellite SAR images, Ground Penetrating Radar data, shallow glacier cores—and the number of land ice bodies analyzed, namely, three tidewater glaciers in Svalbard and one ice cap in Iceland. The main aim is to assess how well popular methods of SAR analysis perform in distinguishing glacier zones, regardless of factors such as the morphologic differences of the ice bodies, or differences in SAR data. We test and validate three methods of glacier zone detection: (1) Gaussian Mixture Model–Expectation Maximization (GMM-EM) clustering of dual-pol backscattering coefficient (sigma0); (2) GMM-EM of quad-pol Pauli decomposition; and (3) quad-pol H/α Wishart segmentation. The main findings are that the unsupervised classification of both sigma0 and Pauli decomposition are promising methods for distinguishing glacier zones. The former performs better at detecting the firn zone on SAR images, and the latter in the superimposed ice zone. Additionally, C-band SAR data perform better than L-band at detecting firn, but the latter can potentially separate crevasses via the classification of sigma0 or Pauli decomposition. H/α Wishart segmentation resulted in inconsistent results across the tested cases and did not detect crevasses on L-band SAR data.

1. Introduction

A growing number of satellites dedicated to Earth observations [1], as well as computational capabilities, offer the possibilities, data and tools to implement high-quality monitoring of the state of the Earth (e.g., [2,3]), including glacier response to climate change. During a time of climate warming, with its strong impact on the Arctic [4,5], glacier monitoring plays an important role in understanding the processes taking place in the cryosphere. In addition, satellite data may be a valuable asset in case of events, such as worldwide lockdowns, when travel restrictions cause gaps in the long time series of glacier measurements [6].

One way of studying a glacier´s condition is by monitoring glacier facies. Glacier facies (e.g., bare ice, superimposed ice, firn) [7,8], also called glacier zones, are parts of a glacier that differ in characteristics such as structure, density, percolation properties or albedo. Together, they form either the accumulation or ablation zone of a glacier (i.e., a zone where a glacier either gains or loses mass in a given time span) [9]. Therefore, changes in the extent of a glacier zone at the end of an ablation season (late summer and autumn) are one of the indicators of glaciers’ state and their response to climate change. Moreover, due to physical differences between glacier zones, information about their extents can support studies of glacier mass balance, hydrology, and other components of the surrounding environment (e.g., [8,10,11]).

Especially valuable for monitoring glacier zones in polar regions are Synthetic Aperture Radar (SAR) satellite data. SAR imagery is insensitive to cloudiness or sunlight. Therefore, unlike multispectral imagery, SAR imagery can be acquired in bad weather conditions and polar night. Additionally, SAR microwaves can penetrate through dry snow, providing information about the glacier’s surface and subsurface [12,13]. As a result, SAR imagery acquired during the accumulation season (winter and spring months) presents the state of a glacier from the end of the ablation season—which is preserved under snow cover due to generally negative air temperatures during the accumulation season.

Several properties of SAR sensors as well as the mode of SAR data acquisition may influence the information provided by SAR imagery, and thus, the results of glacier zone detection. For example, the penetration level of the glacier surface and subsurface depends on the wavelength of the SAR microwave, where longer waves (e.g., L-band SAR sensors) penetrate deeper than shorter waves (e.g., C-band SAR sensors) [14]. Another factor which may lead to differences in the information provided by SAR data is the polarization mode. SAR polarization mode usually varies from one (single-pol), via two (dual-pol), to four SAR channels (quad-pol). With the higher level of polarization mode (e.g., quad-pol), more information on the characteristics of the medium can be retrieved [15,16] and more advanced methods of SAR processing can be applied than in the case of, for example, single-pol SAR imageries. SAR imagery can also be acquired from ascending or descending orbits and under different incidence angles. As a result, the areas of interest can be illuminated differently by SAR sensors if, for instance, the orbit, the incidence angle or the wavelength, differs for the acquisitions.

Information about SAR microwave scattering can be described by different means. One of the most popular coefficients in the analysis of SAR data is the backscattering coefficient (sigma0), which corresponds to the microwave reflectance from the scatterers [17]. In addition, a scattering mechanism can be characterized by using, for example, Pauli or H/A/α decompositions. Pauli decomposition describes contributions of odd-bounce, even-bounce and volume-scattering mechanisms [18], whereas H/A/α represents contributions of entropy, anisotropy and mean alpha angle [19]. Based on the H/A/α, Cloude and Pottier [19] proposed an H/α segmentation—a method of unsupervised classification of the H–α plane into zones that represent different scattering mechanisms (e.g., High Entropy Multiple Scattering or Low Entropy Surface Scatter). Lee et al. [20] improved the method of H/α segmentation by combining a maximum likelihood classifier to the H–α Wishart plane segmentation.

Some successful studies on glacier zone detection and monitoring in polar regions have been carried out, using analysis of sigma0 and comparison to terrestrial data (e.g., [21,22,23,24,25]). These studies are limited to only a few glaciers and to a single- or dual-pol SAR mode of C-band SAR data. Barzycka et al. [25], in a study of changes in glacier zones on Hornsund glaciers (Spitsbergen, Svalbard), recognized limitations in the K-means unsupervised classification of sigma0 of both single- and dual-polarization SAR images in distinguishing glacier zones. The limitations are related to the thickness of glacier facies, heterogeneity of the ice body, and the local character of a structure. Although the classification of glacier zones in that study yielded generally good results, the limitations influenced the classification of, for instance, the superimposed ice (SI) zone in a few of the analyzed cases. On the other hand, a few studies have focused on distinguishing glacier zones using quad-pol SAR data and PolSAR methods, with the assessment based on terrestrial data. For example, Parrella et al. [26,27], and Sharma et al. [28], presented new scattering models validated with airborne PolSAR data from either K-transect (West Greenland) or Austfonna (Nordaustlandet, Svalbard). Doulgeris et al. [29] distinguished glacier zones in Holtedahlsfonna (Spitsbergen, Svalbard) using L-band SAR data from single-, dual- and quad-pol mode and using Wishart and K-Wishart classifiers. Błaszczyk [22] analyzed the possibility of distinguishing glacier zones in Vestfonna (Nordaustlandet, Svalbard) by using both dual-pol C-band SAR data and quad-pol L-band SAR data. The methods used for distinguishing or describing glacier facies in the Błaszczyk [22] study included K-means classification, polarimetric Pauli decomposition, H/A/α and H/α Wishart segmentation. Barzycka et al. [30] expanded the study of Błaszczyk [22] via additional analysis of the backscattering coefficient and a novel application of the objective method of Internal Reflection Energy (IRE) for Ground Penetrating Radar (GPR) data analysis.

The European Space Agency provides free and open SAR data acquired by dual-pol C-band Sentinel-1 satellites with a minimum 6 day revisit time. However, more polarimetric information can be provided by commercial and licensed quad-pol SAR, such as L-band ALOS-2 PALSAR or C-band RADARSAT-2. In commercial solutions, SAR acquisition can be scheduled and the commercial SAR data can have smaller pixel spacing and higher resolution than Sentinel-1. Although commercial SAR data are paid, one high-quality SAR imagery per accumulation season can be sufficient to apply a long-term monitoring of glacier zones [25]. Nevertheless, in this study, we would like to examine if quad-pol SAR data outperforms dual-pol SAR data, or if the results are comparable. If the latter, then the analysis of dual-pol SAR data may be sufficient—and cheaper—for studies related to glacier zones detection.

In this study, we investigate the potential of PolSAR methods for distinguishing glacier zones. We analyze four different land ice bodies, using four satellite quad-pol SAR imageries of either C-band or L-band. The methodology of Błaszczyk [22] and Barzycka et al. [30] is further developed and tested: for quad-pol SAR images, both Wishart H/α segmentation and Pauli decomposition are used. In addition, the results of Pauli decomposition are classified by an unsupervised method of Gaussian Mixture Model–Expectation Maximization (GMM-EM). In order to examine if quad-pol SAR data and PolSAR methods outperform the more popular dual-pol SAR data and sigma0 analysis, similarly to Callegari et al. [31], the dual-pol SAR data are mimicked from quad-pol SAR. The GMM-EM algorithm is used for sigma0 clustering. The quality of the SAR results is assessed based on terrestrial data, namely, GPR measurements and shallow glacier cores. In the case of GPR measurements, two interpretation methods are applied, whose results are compared; namely, a visual interpretation of GPR profiles provided by an expert, and an objective method of unsupervised classification of the Internal Reflection Power (IRP) coefficient.

This study of the effectiveness of PolSAR methods in distinguishing glacier zones is unique in terms of the dataset compiled (four quad-pol satellite SAR images, GPR measurements, shallow glacier cores) and the number of land ice bodies analyzed; namely, three tidewater glaciers located in Hornsund fjord basin (South Spitsbergen, Svalbard) and a land-based ice cap, Langjökull (Iceland). The dataset, as well as the study sites, allows us to assess the potential of the three different methods for distinguishing glacier zones in previous-summer surfaces. We analyze both PolSAR (Pauli decomposition, H/α Wishart segmentation) and backscattering coefficient methods to determine if less accessible quad-pol SAR data outperforms dual-pol SAR in distinguishing glacier zones. We also discuss the potential influence of polarization mode and terrain illumination by SAR on glacier zone detection. Having analyzed glaciers with different characteristics, the results obtained seem to be more generally applicable than those of previous studies of glaciers with similar morphology or climate conditions.

2. Study Sites

In this study, we analyze the glacier zones of three polythermal tidewater glaciers and one temperate land-based ice cap. The tidewater glaciers are Hansbreen, Storbreen and Hornbreen, located in the Hornsund fjord basin (South Spitsbergen, Svalbard), whereas Langjökull, Iceland’s second largest ice cap, is situated in the central western part of the island (Figure 1). Due to the lack of full coverage of Hornbreen by the analyzed SAR image, the area of interest of this glacier is an intersection of SAR area coverage and an ice body outline. In the case of Langjökull, the area of interest is limited both by the footprint area of the SAR imagery and by computational limitations; analysis is therefore provided only for the western part of the ice cap, where terrestrial data were also obtained.

The Icelandic ice cap differs from the Svalbard glaciers in its morphologic form and in its location, being at a lower latitude (Figure 1;

Table 1. Main characteristics of analyzed land ice bodies.

Ice Body	Type	Slope Inclination [°]	Area [km2]	ELA [m a.s.l.]	Velocity [m a−1]
Hansbreen	tidewater glacier	1.7 1	49.4 2	342 3	177 4
Storbreen	tidewater glacier	1.3 1	188.6 2	383 3	132 4
Hornbreen	tidewater glacier	1.3 1	169.4 2	398 3	287 4
Langjökull	ice cap	3.4 5	835 6	1000 5	~75 7

¹ During 2010 [32]. ² During 2017 [25]. ³ During 2014 [36]. ⁴ Average frontal velocity, during 2014 [37]. ⁵ During 2007 [38]. ⁶ During 2019 [39]. ⁷ Maximum horizontal velocity, during 2012 [40].

), affecting the angle of sunshine radiation and the length of insolation. Moreover, as Hornsund and Iceland are located in a polar and subpolar zone, respectively, their climate conditions are different. The mean long-term (1979–2018) annual air temperature in Hornsund is −3.7 °C [34], whilst for the west highlands of Iceland, where Langjökull is located, it is −0.9 °C for the time period 1971−2000 (0.3 °C for 2004−2019) [35]. Mean long-term (1979–2018) annual precipitation in the Hornsund region is 477 mm [34], whereas in the vicinity of Langjökull, it is 727 mm for the period 1971–2000 (724 mm for 2004–2019) [35]. The above-described differences between Hornsund and Langjökull differentially shape their glacier zone state and Equilibrium Line Altitude (ELA [9]; Table 1). Analyzing land ice masses whose morphologic and climatic conditions differ makes this study more generally applicable than studies of glaciers of a single region.

3. Materials and Methods

3.1. Materials

In this section, we describe the satellite and terrestrial data used for distinguishing the glacier zones of Hansbreen, Storbreen, Hornbreen and Langjökull. This includes SAR satellite data (Section 3.1.1), GPR data (Section 3.1.2), and shallow glacier cores (Section 3.1.3).

3.1.1. SAR Data

The presence of water in the snowpack (after e.g., a rain-on-snow event) reduces the SAR penetration depth and influences the scattering mechanism [41,42]. Therefore, we analyzed the quad-pol SAR images acquired by either Advanced Land Observing Satellite-2 Phased Array type L-band Synthetic Aperture Radar-2 (ALOS-2 PALSAR-2) or RADARSAT-2 only during dry snow conditions (

Table 2. Details of SAR images analyzed in this study. Abbreviations: A2—ALOS-2 PALSAR-2 (L-band); SHS QP—Stripmap high-sensitive, Quad Polarization; RS2—RADARSAT-2 (C-band); FQP—Fine Quad Polarization.

Previous Summer Surface	Date of Acquisition	Glacier	Mission; Acquisition Mode	Near Incidence Angle [°]	Far Incidence Angle [°]	Last Positive Temp. Day	Reference Name
2016	10 April 2017	Hansbreen	A2; SHS QP	32.37	35.34	15 March 2017	2016_Hansbreen_A2
2017	12 March 2018	Langjökull	RS2; FQP	38.37	39.84	28 February 2018	2017_Langjökull_RS2
2017		Hansbreen	RS2; Wide FQP	31.72	34.71	28 February 2018	2017_Hansbreen_RS2
	12 March 2018	Storbreen					2017_Storbreen_RS2
		Hornbreen					2017_Hornbreen_RS2
2017	3 April 2018	Hansbreen	A2; SHS QP	17.37	21.90	17 March 2018	2017_Hansbreen_A2

). Assessment of the snow conditions was based on either an on-site inspection or the daily maximum air temperature from the closest weather stations, namely, Hveravellir in Iceland or the Polish Polar Station in Hornsund [43]. As a result, two SAR images of Langjökull (16 March 2018—RADARSAT-2; 17 March 2018—ALOS-2 PALSAR-2) are excluded from the analysis due to a rain-on-snow event shortly before the SAR image acquisition. All SAR data were acquired with ascending orbit and are of the Single Look Complex (SLC) type.

As SAR imagery during the cold season generally presents glacier zone extents from the end of the previous ablation season, in this paper, we refer to the SAR data using the year of the previous summer surface, not the year of acquisition. In addition, we provide the glacier name and SAR mission of the source SAR image; for example, 2017_Langjökull_RS2 is a reference to SAR data for Langjökull’s 2017 summer surface, acquired by RADARSAT-2 in spring 2018. Reference names, as well as details of available SAR quad-pol data, are presented in Table 2.

3.1.2. GPR Data

All GPR data used in this study (

Table 3. Details of GPR measurements used in this study.

Previous Summer Surface	Glacier	Date	Total Length [km]	Sampling Frequency [MHz]	Stacks	Average Distance between Traces [m]
2016	Hansbreen	22 April 2017	100.2	12,791.6	8	1.7
2017	Langjökull	13, 14 March 2018	58.6	12,763.5	4	1.1
2017	Hansbreen	18 April 2018	104.8	16,410.2	4	1.2
2017	Storbreen	26 April 2018	19.7	12,763.5	4	1.2
2017	Hornbreen	26 April 2018	22.8	16,410.2	2	1.1

) were collected with the 800 MHz shielded antenna during dry snow conditions, which were assessed as described in Section 3.1.1. Dry snow conditions were also present between the SAR and GPR data acquisitions. The measurements were collected with the GPR antennas fixed to a sledge, which was pulled behind a snowmobile at a speed of ~20 km h⁻¹. The total length of the GPR profiles and the coverage of the land ice masses by GPR measurements varied depending on safety considerations, the accessibility of the glacier by snowmobiles, and the size of the glacier. Nevertheless, even under the challenging logistics of the GPR surveys, GPR profiles were collected that covered both the accumulation and ablation zones, and that crossed potential glacier facies. The processing of GPR data is described in detail in Barzycka et al. [25].

3.1.3. Shallow Glacier Cores

Shallow glacier cores were collected for snow depth as well as visual analysis of the structure and properties of the glacier surface. The coring was performed on 18 April 2017 (Hansbreen), 22 April 2018 (Hansbreen), 14 March 2018 (Langjökull) and 26 April 2018 (Storbreen and Hornbreen) using a 9 cm diameter corer. Dry snow conditions were preserved between the SAR data acquisition and glacier cores collection. Each of the cores consisted of seasonal snow cover and an underlying structure representing the previous summer’s surface (e.g., ice, SI, firn). The core lengths varied due to differences in snow depth but, in all cases, were sufficient for interpreting the underlying structure. The interpretation was based on the size of the grains, textures, and patterns. Similarly to the SAR and GPR data, we refer to the glacier cores data using the year of the previous summer’s surface, not the year of their collection.

3.2. Methods

Three methods for distinguishing glacier zones, based on SAR satellite data, are studied and compared to the glacier cores and the results of the GPR data analysis (Figure 2). GPR data analysis includes GPR visual interpretation (GPR VI), which is validated by the results of the IRP natural breaks classification. The SAR data analysis methods used in this study are GMM-EM classification of HH/HV sigma0 values, GMM-EM classification of quad-pol SAR, Pauli decomposition, and the H/α Wishart segmentation of quad-pol SAR data. The methods are described, in detail, in Section 3.2.1, Section 3.2.2, Section 3.2.3 and Section 3.2.4.

3.2.1. GPR VI and IRP Natural Breaks Classification

GPR VI was used to identify the glacier zones for further accuracy assessment of the SAR analysis results (Figure 2). The snow cover was identified [44,45] and excluded both from the GPR VI and the calculation of the IRP coefficient. The glacier structure underlying the snow cover was identified based on the reflection characteristics, and classified into the representation of three glacier zones, namely, ice, SI and firn. In this paper, the classes from the terrestrial data analysis (i.e., GPR and shallow glacier cores) are denoted as ice_g, SI_g and firn_g, and represent ice, SI and firn, respectively. The methodology and the results of GPR VI for most of the GPR measurements used in this study are discussed in Barzycka et al. [25].

As a supportive method for GPR VI, the unsupervised natural breaks classification of an IRP coefficient was used to assess the correctness of GPR VI (Figure 2). IRP is an arithmetic mean of reflected energy for a sample within a defined GPR trace time window [46]. This coefficient has been used in glaciological studies of, for example, liquid water presence in glacier systems (e.g., [47,48,49,50]), whereas Barzycka et al. [25] successfully applied the unsupervised classification of IRP in distinguishing glacier zones. Previous studies applying similar coefficients, i.e., the GPR backscatter coefficient [21,51] and Internal Reflection Energy [24,30], to distinguish glacier zones, also yielded good results.

Barzycka et al. [25] recommended the natural breaks classification of the IRP as an alternative to, or a supportive method for, GPR visual interpretation. As the methodology and the results of IRP natural breaks classification, for most of the GPR measurements used in this study, are discussed at length in Barzycka et al. [25], here, we present the accuracy assessment of IRP natural breaks classification as a measure of agreement with the GPR VI results. The accuracy assessment consists of precision (user’s accuracy), recall (producer’s accuracy), F-score, and Kappa metrics [52]. The high values of those scores indicate that the objective IPR natural breaks classification supports the results of GPR VI.

3.2.2. GMM-EM Classification of Dual-Pol SAR Sigma0

Sigma0 is a quantity that corresponds to the microwave reflectance from the scatterers that was recorded by the SAR sensor [17]. In studies by Błaszczyk [22] and Barzycka et al. [30], one of the tested methods for distinguishing glacier zones is based on the sigma0 coefficient of dual-pol SAR data. In this study, however, we assess the possible underperformance of the sigma0 dual-pol SAR data compared to the quad-pol PolSAR methods. To enable such a comparison, sigma0 is retrieved from two channels of quad-pol SAR data [31]. Therefore, unlike Błaszczyk [22] and Barzycka et al. [30], the unsupervised classification of dual-pol SAR sigma0 should be comparable to quad-pol SAR methods, as all sets of SAR data have common orbits or incidence angles of SAR acquisition. The polarization channels retrieved from quad-pol SAR data are HH and HV. This choice is based on the dual-pol SAR polarization modes that were available for the Hornsund area in the Barzycka et al. [25] study.

The algorithm of the unsupervised classification of sigma0, which is used in this study, is GMM-EM [53]. GMM-EM has been successfully applied in many different studies (e.g., [54,55,56,57]). In this study, GMM-EM is chosen as an alternative to a popular K-means unsupervised classification, which was also used in the studies of Błaszczyk [22] and Barzycka et al. [25,30]. GMM-EM was chosen in preference to K-means because the decision boundaries of the GMM-EM clustering depend on the mean and variance of the data, whereas in K-means, they are shaped based on the distance metric from the cluster’s centroids [53]. Therefore, GMM-EM is more flexible than K-means and can be applied to clusters of more complex patterns in data [58], and that are, possibly, also areas of heterogeneity of ice or SI of a local character, which were not detected by the K-means method in [25].

As a pre-processing step for sigma0 retrieval (Figure 2), the SAR data was calibrated [59,60] and filtered using a Refined Lee filter with a 11 × 11 window size, a 3 × 3 target window size, and a single number of looks [18]. A terrain correction algorithm was then applied [59,61] using ArcticDEM-based Digital Elevation Models [33,62]. In order to mimic dual-pol SAR satellite data, only sigma0 values for HH and HV polarization were chosen and scaled to decibels. All SAR pre-processing steps were performed in SNAP 8.0 software [63].

As an input for the GMM-EM clustering, only sigma0 values within glacier outlines (Figure 1) were considered. The algorithm used for the clustering was GMM-EM, with a general covariance matrix for each component, 50 initializations, and 1000 maximum iterations with a default early stopping of a 0.001 convergence threshold [64]. In order to avoid subjectivism, a number of clusters (components) for each dataset were chosen based on Silhouette [65], Davies–Bouldin [66] and Calinski–Harabasz [67] scores, as well as clustergams [68,69]. For computational reasons, up to 10 components were scored by the above algorithms. Similarly to Barzycka et al. [25], the components were later carefully grouped to represent glacier zones (referred to as ice_SAR, SI_SAR, firn_SAR) based on terrestrial data, UAV flights, and satellite optical imageries.

For the accuracy assessment of the sigma0 method, categorical values of ice_SAR, SI_SAR and firn_SAR were retrieved and compared to the results of GPR VI. Metrics used for the accuracy assessment were precision, recall and F-score. The resulting scores are described in detail, as even a few hundredths of difference may impact the detection of glacier zone boundaries and the assessment of, for instance, the glacier zone areas. In the rest of the text, we refer to the GMM-EM classification of dual-pol SAR sigma0 also as sigma0+GMM-EM.

The Geopandas 0.8.1 Python package [70] was used to retrieve input data for the clustering within glacier outlines. The GMM-EM, metrics calculation, and scoring were performed using the ScikitLearn 0.23 [64] and Clustergram 0.5.0 [68,69] Python packages.

3.2.3. GMM-EM Classification of Quad-Pol SAR Pauli Decomposition

Pauli decomposition [18] of quad-polarimetric SAR data is one of the most popular methods for analyzing scattering mechanisms, also used in cryospheric studies (e.g., [22,30,71,72,73,74]). It provides information on the contribution of odd-, even-bounce and volume-scattering mechanisms for each pixel, which can be related to different glacier zones. Błaszczyk [22] and Barzycka et al. [30] are examples of the successful application of Pauli decomposition in the description of glacier zones compared to GPR VI. In those studies, however, Pauli decomposition was not classified in any manner (supervised or unsupervised). Therefore, clear extents of the glacier zones were not presented, and an accuracy assessment based on terrestrial data was not performed. In this study, the results of Pauli decomposition are classified in an unsupervised way by GMM-EM (see Section 3.2.2.) to retrieve the extents of the glacier zones and perform an accuracy assessment based on GPR VI. GMM-EM classification of quad-pol SAR Pauli decomposition is referred to here as Pauli+GMM-EM.

As a pre-processing step for Pauli decomposition, all SAR images were calibrated [59,60] (Figure 2). A T3 polarimetric matrix was then generated [18] and filtered using the Refined Lee filter [18], as described in Section 3.2.2. Based on this, Pauli decomposition was performed [18] and the data were geometrically corrected [59,61] using the DEMs described in Section 3.2.2. Results were scaled to decibels. All described steps were performed in SNAP 8.0 software [63]. In the next step, the pixels for all Pauli decompositions were retrieved within each glacier outline (Figure 1) and classified into firn_SAR, SI_SAR or ice_SAR, using GMM-EM in the same manner as described in Section 3.2.2. Results were assessed based on the GPR VI results using precision, recall and F-score metrics.

3.2.4. H/α Wishart Segmentation of Quad-Pol SAR Data

H/α Wishart segmentation is a method of unsupervised classification based on the H–α plane and the Complex Wishart Classifier [19,20,75]. The H-α plane is divided into nine zones, eight of which are associated with a different scattering mechanism (one is a non-feasible region) [19]. Due to fixed boundaries of the H–α plane and the locations of the final class centers therein, this method is particularly suitable for application as an unsupervised method for distinguishing glacier zones, where each glacier zone could potentially be represented by one, or a few, particular H–α scattering mechanisms. Studies such as Błaszczyk [22] and Barzycka et al. [30] show that glacier zones can be distinguished with H–α Wishart segmentation. However, there is no previous study that applies this method for several SAR images or glaciers and assesses the accuracy of the results based on terrestrial data. In this study, we applied H/α Wishart segmentation to distinguish the glacier zones of four land ice masses and assessed its performance based on terrestrial data.

To achieve H/α Wishart segmentation, each of the SAR data was multilooked (number 1 and 2 of the range and azimuth looks, respectively) and a T3 polarimetric matrix was generated [18] (Figure 2). Next, based on the T3 matrix, H/α decomposition [75] was performed, followed by H/α Wishart segmentation [20] with 100 maximum iterations and a 5 × 5 window size for the Boxcar filter parameter. Due to software limitations, the H/α Wishart segmentation was applied to a full SAR imagery, not an AOI subset (Figure 1). As a final step, a terrain correction was applied [59,61]. The described process was performed in PolSARPro 6.0 software [76], and the workflow is based on [77]. The results of the H/α Wishart segmentation were grouped into SAR glacier zones, and the accuracy assessment was performed as described in Section 3.2.2.

4. Results

The results of this study are focused on the analysis of terrestrial data (Section 4.1) and on distinguishing glacier zones based on SAR data (Section 4.2). The outcomes of the tested SAR methods, i.e., the GMM-EM classification of sigma0 and of Pauli decomposition as well as H/α Wishart segmentation, are reported separately in Section 4.2.1, Section 4.2.2 and Section 4.2.3.

4.1. Analysis of Terrestrial Data

Results of the GPR VI show that the profile analyses are spatially consistent in all study sites (Figure 3 and Figure 4). Ice_g and firn_g are generally identified in the glaciers’ lower and upper altitudes, respectively. SI_g—if present—is identified either between the ice_g and firn_g classes or in short sections between the ice_g. The above is in good agreement with what is known about glacier facies present in the glacier system [7,8], whereas short sections of SI_g between ice_g are formed locally due to glacier topography [25].

The metrics of IRP natural breaks classification for ice_IRP and firn_IRP are generally high (

Table 4. Results of the accuracy assessment of the IRP natural breaks classification [25].

Previous Summer Surface	Glacier	Class	Precision	Recall	F-Score	Kappa
2016	Hansbreen	iceIRP	1.00	0.95	0.97	0.96
		firnIRP	0.85	0.99	0.91
2017	Langjökull	iceIRP	0.99	0.93	0.96	0.96
		firnIRP	0.94	0.99	0.97
		iceIRP	0.98	0.98	0.98
2017	Hansbreen	SIIRP	-	0.00	-	0.96
		firnIRP	0.92	0.99	0.95
		iceIRP	0.98	0.94	0.96
2017	Storbreen	SIIRP	0.89	0.98	0.93	0.96
		firnIRP	1.00	0.98	0.99
		iceIRP	0.99	0.99	0.99
2017	Hornbreen	SIIRP	-	0.00	-	0.96
		firnIRP	0.86	1.00	0.92

). The exception is the precision score for firn_IRP in both the Hansbreen 2016 summer surface (0.85) and the Hornbreen 2017 summer surface (0.86). The former is mainly due to a strong noise in the GPR data collected with a low battery voltage, which influenced the natural breaks thresholds [25]. The latter is due to the misclassification of a shallow firn_g as ice_IRP, as shallow firn_g is characterized by lower values of IRP [25]. SI_IRP was classified only in the case of the Storbreen 2017 summer surface, with generally good results. The lack of representation of SI_g in the IRP results for Hansbreen and Hornbreen for the 2017 summer surface is most likely due to the local character of this glacier’s structure [25]. Thus, the number of samples of SI_g was not enough to distinguish it as a separate SI_IRP class but, at the same time, it did not influence the overall accuracy of the IRP classification (Kappa score of 0.96 for both Hornbreen 2017 and Hansbreen 2017). Overall, the good results of the IRP accuracy assessment show a strong agreement of the subjective method of GPR VI with the objective method of IRP, which at the same time, indicates a high quality of GPR VI.

In general, the information retrieved from the shallow glacier cores was in good agreement with the results of the GPR VI (Figure 3 and Figure 4). Therefore, ice_g and firn_g were identified in shallow glacier cores in the glaciers’ lower and upper parts, respectively. Glacier cores retrieved close to the firn line of Hansbreen for the 2016 and 2017 summer surface consisted of polycrystalline glacier ice with firn inclusions [25]. In the same area, a relatively high number of scattering elements were identified in the GPR profiles. Due to the predominance of glacier ice, this was classified as ice_g; however, the occurrence of incompletely developed polycrystalline glacier ice is noted and referred to as a “transition area” (or as ice + firn_g). Two shallow glacier cores in Hornbreen are in disagreement with the GPR VI of the nearest profiles. This, however, can be explained by local depressions in the glacier topography, and a possible influence of positioning error for a handheld GNSS device used during glacier core drilling [25].

4.2. Distinguishing Glacier Zones Based on SAR Data

In Section 4.2.1, Section 4.2.2 and Section 4.2.3, we describe, in detail, the results of distinguishing glacier zones based on SAR data using the GMM-EM classification of sigma0, GMM-EM classification of Pauli decomposition, and the H/α Wishart segmentation methods. We refer to the spatial characteristics of the distinguished zones (Figure 3 and Figure 4) as well as to scores for recall, precision, and the F-score metrics (Figure 5) for each analyzed ice body and method.

4.2.1. GMM-EM Classification of Dual-Pol SAR Sigma0

GMM-EM clustering of sigma0 based on HH and HV generally yields good results for glacier zone detection on Hansbreen (Figure 3). In the case of the detection of firn_g at this glacier in the ALOS-2 PALSAR-2 imageries (2016_Hansbreen_A2, 2017_Hansbreen_A2), misclassification occurs especially in the area of firn_SAR and ice_SAR boundaries, where shallow firn occurs and this is classified as ice_SAR. Such a pattern is less visible in the 2017_Hansbreen_RS2 case. This is reflected in the recall values for firn_SAR, where 2016_Hansbreen_A2 and 2017_Hansbreen_A2 received 0.82 and 0.84 recall scores, respectively, whereas 2017_Hansbreen_RS2 was scored at a 0.88 level of recall (Figure 5). In addition, and in contrast to the ALOS-2 PALSAR-2 results, in 2017_Hansbreen_RS2, the transition area is represented by firn_SAR patches below the firn line. Firn_g in Tuvbreen (Hansbreen’s tributary glacier, Figure 1) is detected in the ALOS-2 PALSAR-2 and RADARSAT-2 imagery. However, the area of this class at Tuvbreen is greater in clustering results of RADARSAT-2 than in ALOS-2 PALSAR-2 imagery; this difference could also impact the classification’s metrics. SI_g, labelled for the 2017 summer surface at Deileggbreen (Hansbreen’s tributary glacier, Figure 1), was not detected by the sigma0+GMM-EM method, and the algorithm classified it as ice_SAR at both 2017_Hansbreen_A2 and 2017_Hansbreen_RS2. In addition, the area of ice_SAR in the upper part of Hansbreen (at ~77°8′ N) is larger for the ALOS-2 PALSAR-2 results than for RADARSAT-2. In addition, sigma0+GMM-EM of ALOS-2 PALSAR-2 imageries made it possible to distinguish a separate class for highly crevassed areas, such as at the front of Hansbreen.

Firn_SAR at Storbreen is in good agreement with firn_g labels (Figure 4). This is reflected in the high recall and precision scores for firn detection (0.96 and 0.99, respectively, Figure 5). The SI_SAR class was distinguished at Storbreen with a relatively high score for precision (0.91) and a moderate result for recall (0.71). The latter is mainly due to the misclassification of SI_g in the area between Glimseen and Drevbreen (tributary glaciers of Storbreen, Figure 1). SI_SAR is also detected in the lower parts of the ablation area in the shape of patches. GPR measurements, however, do not cover those areas; therefore, it has neither positive nor negative influence on the scores and the final assessment of the method. The ice_SAR class for Storbreen has one of the lowest scores for precision (0.84), which may be a result of the classification of SI_g as ice_SAR in the area between Glimseen and Drevbreen. The sigma0+GMM-EM of 2017_Storbreen_RS2 does not result in a separate class for crevasses.

The boundary of firn_SAR for the Hornbreen glacier (Figure 4) is located higher than in the case of firn_g (recall score at 0.76 level, Figure 5). Similarly to the Hansbreen cases, a relatively small class of SI_g is not detected for Hornbreen by the sigma0+GMM-EM method. The precision score of 0.90 for ice_SAR at Hornbreen is most likely due to the misclassification of the glacier zones (firn_g, SI_g) in favor of ice_SAR in the vicinity of the firn_SAR–ice_SAR boundary. The clustering of 2017_Hornbreen_RS2 does not result in a separate class, which could represent the glacier’s crevasses; instead, there is a relatively large area of frontal crevasses classified as firn_SAR. The misclassification of this heterogeneity of the ice body does not influence the classification scores, as this area was not covered by GPR measurements, for safety reasons.

In the case of Langjökull, the firn_g class is misclassified in the vicinity of the firn_SAR line in favor of the ice_SAR class (Figure 4). This pattern is especially visible on the southern slopes of the ice cap, where the firn_SAR boundary is higher than the firn_g labels. On the northwest slope, such a pattern is not present: firn_SAR and firn_g are in good agreement. Furthermore, a small patch of firn_g (~64°38′N, 20°30′W) is correctly classified as firn_SAR on the northwest slope of the ice cap. There is, however, misclassification of ice_g on the northwest slope as firn_SAR in the ablation area (~ 64°38′–64°40′N, 20°30′–20°33′W). In addition, singular pixels of ice_SAR within the main firn_SAR area are dominant on the slopes facing south. As a result, the overall scores for precision and recall for firn_SAR at Langjökull are 0.97 and 0.84, respectively, whereas for ice_SAR, they are 0.83 and 0.97, respectively (Figure 5). The firn_SAR area, which does not influence the metrics (no terrain data collected) but is worth noting, is located south of Flosajökull (outlet glacier of Langjökull, Figure 1). This firn_SAR, which is discontinuous from the firn_SAR covering the ice cap dome, is formed due to local depressions as well as the presence of narrow steep valley outlet glaciers. Sigma0+GMM-EM for Langjökull does not result in a separate class that could represent crevasses.

4.2.2. GMM-EM Classification of Quad-Pol SAR Pauli Decomposition

The firn_SAR line at both 2016_Hansbreen_A2 and 2017_Hansbreen_A2 (Figure 3) is detected higher than indicated by firn_g (recall score is 0.73 for both cases of firn detection, Figure 5) using the Pauli+GMM-EM method, whereas in the case of 2017_Hansbreen_RS2, the firn_SAR line is in good agreement with the terrestrial data results (recall score: 0.82) and a few patches of firn_SAR below the firn_SAR line are also distinguished, indicating glacier ice that is not fully developed (transition area). Firn_SAR at Tuvbreen is detected on 2016_Hansbreen_A2, 2017_Hansbreen_A2 and 2017_Hansbreen_RS2 with good agreement with firn_g. A rather local patch of SI_g at Deileggbreen is not distinguished as a separate class in either 2017_Hansbreen_A2 or 2017_Hansbreen_RS2. Similarly to the sigma0+GMM-EM results, the ice_g patch in the upper part of Hansbreen is distinguished in the Pauli+GMM-EM results of both the ALOS-2 PALSAR-2 and RADARSAT-2 images; the patch, however, covers a larger area in the ALOS-2 PALSAR-2 than RADARSAT-2 results. By applying Pauli+GMM-EM onto the ALOS-2 PALSAR-2 images, it is possible to distinguish a separate class of highly crevassed areas.

The extent of firn_SAR in the results of Pauli+GMM-EM for 2017_Storbreen_RS2 (Figure 4) is in good agreement with firn_g (recall score: 0.97, Figure 5). SI_SAR is also detected in this glacier, including the area between Glimseen and Drevbreen. Therefore, the SI_SAR recall score (0.86) of Pauli decomposition and clustering is higher than for the sigma0+GMM-EM method by 0.15. Similarly to the sigma0+GMM-EM, this method also distinguished SI_SAR patches in the ablation area, which are difficult to assess due to the lack of terrestrial data in those parts of the glacier. In the case of ice_SAR, the precision and recall scores for 2017_Storbreen_RS2 are relatively high (0.93 and 0.98, respectively). However, there is no separate class for crevasse representation.

Similar to the sigma0+GMM-EM results, the boundary between firn_SAR and ice_SAR at 2017_Hornbreen_RS2, classified by the Pauli+GMM-EM method, is higher than indicated by the terrestrial data (Figure 4). As a result, the recall score of firn_SAR is at a 0.72 level (Figure 5). In addition, no SI_SAR is detected, and the small area covered by SI_g is classified as ice_SAR (therefore, the precision score of ice_SAR is 0.90). The Pauli+GMM-EM method of the RADARSAT-2 image for the Hornbreen 2017 summer surface did not allow a separate class to be distinguished for vast areas of frontal crevasses at this tidewater glacier. Instead, the area of Hornbreen that is of heterogeneous character was classified as firn_SAR.

The firn_SAR line at Langjökull’s southern slopes is detected higher in the Pauli+GMM-EM results than indicated by the terrestrial data (Figure 4). Additionally, in the firn_SAR class, singular pixels of ice_SAR are present, mainly on the southern slopes. In the northwest part of the ice cap, the extent of firn_SAR is in good agreement with firn_g; however, a small patch of firn_g (~ 64°38′N, 20°30′W) is not detected. The firn_g area between Langjökull dome and Geitlandsjökull (lateral glacier of Langjökull, Figure 1) is classified as a mix of firn_SAR and ice_SAR. Overall, due to the misclassification of firn_g as ice_SAR, the recall score for firn_SAR and the precision score for ice_SAR for 2017_Langjökull_RS2 are both at a level of 0.79. The score of one for both the precision of firn_SAR and the recall of ice_SAR indicates that the misclassification of ice_g as firn_SAR, in this case, is rare. This is supported by the lack of firn_SAR in the northwest ablation area (~ 64°38′–64°40′N, 20°30′–20°33′W), in contrast to the sigma0+GMM-EM results, for example. The firn_SAR area in the south of Flosajökull is not present.

4.2.3. H/α Wishart Segmentation of Quad-Pol SAR Data

The firn_SAR line of 2016_Hansbreen_A2 is detected lower when using the H/α Wishart method than indicated by the terrestrial data (Figure 3). Additionally, large areas in the ablation area where crevasses are present are classified as firn_SAR, not ice_SAR. The ice_g patch in the upper part of Hansbreen is detected in the form of a few pixels, not a patch. All of this produces a 0.89 precision score for firn_SAR of 2016_Hansbreen_A2 (Figure 5). Contrary to the results of 2016_Hansbreen_A2, the firn_SAR line of 2017_Hansbreen_A2 is detected higher than indicated by the terrestrial data. Several distinct ice_SAR patches are visible in the accumulation area, and crevassed areas close to the front of the glacier are classified as firn_SAR. Those misclassifications are reflected in the clustering scores, where firn_SAR of 2017_Hansbreen_A2 with the H/α Wishart method had a 0.73 score for precision for firn_SAR and 0.83 for recall, whereas ice_SAR scored 0.93 for precision and 0.91 for recall. Neither the SI_g patch at Deileggbreen nor highly crevassed areas are classified as a separate class in 2017_Hansbreen_A2. In the case of the results of H/α Wishart for 2017_Hansbreen_RS2, the firn_SAR line is generally detected in good agreement with the terrestrial data. In addition, both the ice_g patch in the upper part of the glacier and the firn_g patch at Tuvbreen are well detected in this case. However, neither the SI_g patch at Deileggbreen nor the frontal crevasses are distinguished as a separate class in 2017_Hansbreen_RS2. Overall, the firn_SAR recall and precision scores of 2017_Hansbreen_RS2 are 0.84 and 0.95, respectively.

In the case of 2017_Storbreen_RS2 classified by the H/α Wishart method, the firn_SAR line is higher than indicated by the terrestrial data (Figure 4), which could influence the firn_SAR recall score (0.94, lower than for sigma0+GMM-EM and Pauli+GMM-EM, Figure 5). SI_SAR is represented by a separate class; however, the recall score of SI_SAR (0.56) indicates that a significant part of this area is misclassified as ice_SAR (precision of ice_SAR: 0.77). Part of this misclassification occurs in the area between Glimseen and Drevbreen. Moreover, the SI_SAR patches in the lower part of the glacier in the ablation area, which are detected by the sigma0+GMM-EM and Pauli+GMM-EM methods, are less distinct in the case of the H/α Wishart classification. Frontal crevasses are not classified as a separate class but as a firn_SAR.

H/α Wishart distinguishes the firn_SAR line of 2017_Hornbreen_RS2 higher than indicated by the terrestrial data and higher than the results of the other analyzed methods (Figure 4). This misclassification of the firn_SAR as ice_SAR in the vicinity of the firn line is reflected in the firn_SAR recall score, 0.69 (Figure 5). Moreover, similarly to other methods, SI_g is not detected at 2017_Hornbreen_RS2. Vast, crevassed areas of Hornbreen’s ablation zone are, similarly to other methods, classified as firn_SAR. However, this misclassified area of firn_SAR at the front of the glacier is smaller compared to the results of either sigma0+GMM-EM or Pauli+GMM-EM for 2017_Hornbreen_RS2.

As in the results of previous methods, the firn line in the southern slopes of 2017_Langjökull_RS2 is located higher than in the northwest part of the ice cap compared to the terrestrial data (Figure 4). The H/α Wishart method also distinguishes a small patch of firn_g ~ 64°38′N, 20°30′W, and there is a distinct firn_SAR between Geitlandsjökull and Langjökull dome. In general, the firn_SAR line is in better agreement with the GPR VI at both the south and northwest slopes than sigma0+GMM-EM or Pauli+GMM-EM; and thus, the firn_SAR recall score for 2017_Langjökull_RS2 is 0.91. Similar to the sigma0+GMM-EM clustering results, in the ablation area of the ice cap (~64°38′–64°40′N, 20°30′–20°33′W), a misclassification of ice_g as firn_SAR occurs. In the area south of Flosajökull, the firn_SAR area is detected.

5. Discussion

This is the first study of the possibility of glacier zone detection based on more than one quad-pol SAR image (including both C-band RADARSAT-2 and L-band ALOS-2 PALSAR-2), with the quality assessment based on terrestrial data. The results of the tested methods are generally similar for both Hornsund glaciers and Langjökull, despite the differences in morphology and climate conditions of the regions. However, the results also show that the performance of each method, or type of SAR data, varied across all of the studied cases. Therefore, in studies focused on the assessment of various data or methodologies (e.g., [30]), it is important to test performance in distinguishing glacier zones on more than one SAR imagery of a particular type. Additionally, good-quality terrestrial data are an important asset in methodological studies and allow performance to be compared.

Because of the deeper penetration of L-band (ALOS-2 PALSAR-2) compared to C-band (RADARSAT-2) in the glacier’s body [14,78], we observed the higher location of the firn_SAR line, or larger area of the ice_SAR patch at ~77°80′N, in both sigma0+GMM-EM and Pauli+GMM-EM results of 2016_Hansbreen_A2 and 2017_Hansbreen_A2 in comparison to the 2017_Hansbreen_RS2 results. This is also reflected in the lower recall scores of firn_SAR for the above-mentioned ALOS-2 PALSAR-2 cases.

The firn_SAR of all tested methods was in better agreement with the terrestrial data results on the northwest slope than the southern slopes of Langjökull. Additionally, singular pixels of ice_SAR within the firn_SAR classes were present, mainly on the southern slopes of Langjökull. This could be due to, e.g., differences in morphological or environmental conditions between the slopes that shape the ice cap structure, an influence of SAR speckle filtering, or a difference in terrain illumination by microwaves during the SAR image acquisition [74,79,80]. More SAR and terrestrial data are needed to define the origin of the differences in glacier zone detection on the slopes of Langjökull. However, in order to reduce the potential influence of the SAR acquisition on how well glacier zones are distinguished in the future, solutions suggested by Callegari et al. [31] could be implemented, i.e., either a classification based on SAR data acquired from both ascending and descending orbits or implementation of the local incidence angle information on non-topographic corrected SAR data classification.

Barzycka et al. [25] found problems in distinguishing SI_g based on dual-pol data. Here, despite more polarimetric information being provided by quad-pol than dual-pol SAR images, the SI_g of Hansbreen and Hornbreen was not recognized as a separate class by any of the methods in this study, or by using different SAR band lengths. This can be due to (1) the thin layer of SI_g [25], which does not have a large influence on overall SAR reflectance in either C- or L-band; (2) low representation, or poor separation, of data representing the local SI_g, so it is not recognized by the GMM-EM algorithm [81]; or (3) too strong a speckle filtering algorithm [82,83], so that the local SI_g is averaged with other zones.

Another issue of SI_g distinction is related to the unknown origin of SI_SAR patches in the low ablation zone of Storbreen. The SI_SAR patches occur in the results of all tested methods. Laska et al. [36] indicate that the low slopes, low number of crevasses, and poor englacial drainage system of Storbreen promoted the formation of meltwater ponds in the 2014 ablation season. This could explain the pattern of SI_SAR patches, whose origin could be a frozen meltwater pond from the ablation season, which was not drained out from the glacial system before the accumulation season. In the Barzycka et al. [25] study, such patches were reclassified as ice_SAR. From the perspective of the monitoring of glacial state or mass balance assessment, the reclassification of such patches can still be a good choice—depending on whether the structure was present throughout the glaciological year. Additional GPR measurements or shallow glacier core drilling in the ablation zone of Storbreen would be beneficial to confirm the potential origin of the SI_SAR patches as frozen meltwater ponds.

Barzycka et al. [25] reported another issue in distinguishing glacier zones, namely, the misclassification of highly crevassed areas of glacier ice as firn_SAR. This was especially problematic in the case of Hornbreen glacier zone detection. This glacier is characterized by a higher frontal velocity than the other analyzed glaciers (Table 1) and the high heterogeneity of glacier ice (caused by a crevassed surface), which—as a result of the K-means classification of sigma0 of Sentinel-1 (C-band) —was misclassified as firn_SAR in Barzycka et al. [25]. In this study, a separate class for crevassed areas, resulting from ALOS-2 PALSAR-2 data and the sigma0+GMM-EM and Pauli+GMM-EM methods, indicates that SAR data from longer-wavelength sensors (e.g., L-band) could potentially resolve the issue of the misclassification of high heterogeneity glacier ice. This, however, should be tested by glacier zone detection on, for example, Hornbreen using L-band SAR data. As the H/α Wishart method did not result in a separate class for the crevassed areas, the method of SAR data processing and classification should be carefully chosen. The possibility of distinguishing, for instance, crevassed areas using L-band SAR data was also noticed by König et al. [84]. Alternatively, supervised methods of crevasse detection can be applied to differentiate between firn_SAR and crevasses [85].

Although there is no clear winner when it comes to comparing performance, regarding distinguishing glacier zones by either the sigma0+GMM-EM, Pauli+GMM-EM or H/α Wishart methods, the results indicate that some of the methods perform better at firn or SI detection than the others. For example, recall scores for firn detection based on Pauli+GMM-EM are generally lower than by other methods, which means that firn_g is generally more often misclassified in Pauli+GMM-EM than in sigma0+GMM-EM, for example. That firn_SAR received generally higher recall scores in sigma0+GMM-EM than Pauli+GMM-EM indicates that the scattering characteristics described by Pauli polarimetric decomposition are less distinct for shallow firn than its backscattering coefficient. On the other hand, SI_g on Storbreen was best represented by the Pauli+GMM-EM method in this study, possibly due to the distinct scattering characteristic of this medium represented by Pauli decomposition, which was then well separated by GMM-EM. However, additional analysis of the different glaciers where SI_SAR was detected would be needed to support recommending this method for SI detection. Finally, the H/α Wishart method gave inconsistent results across the tested cases, so it is difficult to assess the performance of this method. For example, the overall scores for the 2016_Hansbreen_A2 case are very good compared to other methods, but this is not the case for 2017_Hansbreen_A2. Moreover, the results of H/α Wishart for both 2016_Hansbreen_A2 and 2017_Hansbreen_A2 (Figure 3) suggest a rather significant change in the firn_SAR zone after one season; this, however, is not supported by the terrestrial data results. The reason for these different results may be (1) the fixed boundaries of the H–α plane in which the cluster centers are located, whereas GMM-EM classifies the data based on means and variations; (2) the applied workflow (Figure 2), where H/α Wishart is not limited to AOIs, so the segmentation is performed on full SAR imagery, not only on glacier extent. Potentially, the different results of H/α Wishart for 2016_Hansbreen_A2 and 2017_Hansbreen_A2 could also be related to the influence of the incidence angle of the SAR data acquisition, as the incidence angle of 2016_Hansbreen_A2 data is in the range 32.37°–35.34°, whereas for 2017_Hansbreen_A2 it is 17.37°–21.90° (Table 2). This would mean that either H/α Wishart is more sensitive to the image acquisition settings than the other tested methods, or that the differences in the overall processing of H/α Wishart (different software and workflow than sigma0+GMM-EM or Pauli+GMM-EM, Figure 2) influenced the results. H/α Wishart is also the only method with which large crevasses are not detected as a separate class in L-band SAR data—this could also be due to either the location of the cluster centers within the fixed boundaries of the H–α plane or the lack of an AOI subset in the applied workflow. Nevertheless, despite either the fixed boundaries of the H–α plane or the differences in data acquisition or processing, the overall results of H/α Wishart are generally satisfactory for glacier zone detection.

The comparable performance of sigma0+GMM-EM to quad-pol Pauli+GMM-EM or H/α Wishart shows that, in the analyzed cases, polarimetric methods—based on less accessible quad-pol SAR data—do not outperform an unsupervised classification of a popular backscattering coefficient. This is especially promising for the application of free and open Sentinel-1 data in distinguishing glacier zones [25]. Nevertheless, the potential of quad-pol SAR satellites in distinguishing glacier zones can be further studied by, e.g., the classification of multiple polarimetric decompositions as a single dataset [74], or by the classification of higher resolution commercial SAR datasets [86] in comparison to lower resolution open SAR data [87]. Based on the results of this study, we recommend that further analyses are performed based on several SAR imageries and with good-quality terrestrial data.

6. Conclusions

This study compared the performance of glacier zone detection methods for four different land ice masses using quad-pol SAR imageries of C- and L-band. The methods tested were (1) Gaussian Mixture Model–Expectation Maximization unsupervised classification of HH/HV sigma0; (2) Gaussian Mixture Model–Expectation Maximization unsupervised classification of quad-pol Pauli decomposition; and (3) quad-pol H/α Wishart SAR segmentation. The land ice masses analyzed were three glaciers of the Hornsund fjord basin in Svalbard (Hansbeeen, Storbreen, Hornbreen) and one ice cap in Iceland (Langjökull), which differ both in morphologic and climatic conditions, making this study more generally applicable than studies focused on the glaciers of a single region. The results of glacier zone distinction by SAR analysis were validated by terrestrial data (shallow glacier cores, Ground Penetrating Radar). The main conclusions are:

• The results of the unsupervised classification (Gaussian Mixture Model–Expectation Maximization algorithm) of both HH/HV sigma0 and Pauli decomposition are the most promising for distinguishing glacier zones.
• Firn on analyzed SAR images is better represented by the classification results of dual-pol sigma0 than by quad-pol Pauli decomposition and classification.
• Better results for the detection of the SI of Storbreen were obtained by the unsupervised classification of quad-pol Pauli decomposition than of dual-pol sigma0. However, to confirm that the unsupervised classification of Pauli decomposition performs better than other methods in distinguishing SI, more tests on different glaciers are needed.
• The H/α Wishart method gave less satisfactory results than the unsupervised classification of either sigma0 or Pauli decomposition. This is due to inconsistent results with regard to distinguishing glacier zones on Hansbreen, which were assessed based on terrestrial data and accuracy metrics. The inconsistency in the H/α Wishart results is probably determined by either the fixed boundaries of the H–α plane where the cluster centers are located or by differences in the processing workflow in comparison to the unsupervised classification of sigma0 or Pauli.
• To detect a firn zone on SAR images, shallower-penetrating C-band RADARSAT-2 data give better results than L-band ALOS-2 PALSAR-2 when the unsupervised classification of either sigma0 or Pauli decomposition is used.
• The unsupervised classification of dual-pol sigma0 is not outperformed by the results of the classification of quad-pol SAR data and polarimetric methods. This is especially promising in terms of the better availability of dual-pol than quad-pol SAR data.
• The heterogeneity of the glacier ice body could potentially be distinguished by L-band SAR data and the application of the unsupervised classification of either sigma0 or Pauli decomposition. To support this, more tests are needed, especially for glaciers with highly crevassed areas.
• Despite the differences in morphology or climate conditions of the land ice masses of Svalbard and Iceland, the assessed quality of the results of the tested methods are comparable.