Retrieval of Spatial and Temporal Variability in Snowpack Depth over Glaciers in Svalbard Using GPR and Spaceborne POLSAR Measurements

Abstract

The highly dynamic nature of snow requires frequent observations to study its various properties. Keeping this in mind, the present investigation presents results from the analysis of fully polarimetric synthetic aperture radar (POLSAR) parameters for the development of a snow depth (SD) inversion model for SD retrieval. Snow depth retrieved using ground penetrating radar (GPR) at 500 MHz over Austre Grønfjordbreen in the Svalbard region was used to understand the behaviour of certain polarimetric parameters. A significant correlation was found between field-measured SD and POLSAR parameters, namely coherence and normalized volume scattering power (R² = 0.84 and R² = 0.73, respectively.) Using the POLSAR scattering powers obtained from the six-component model-based decomposition (6SD), the heterogeneity and anisotropic behaviour in the firn areas are also explained. Further, based on the analyses shown in this work, a polarimetric parameter-based SD inversion algorithm have been proposed and validated. The univariate model with co-polarization coherence has the highest correlation (R² = 0.84, Root Mean Square Error (RMSE) = 0.18). We have even tested several multivariate models for the same, to conclude that a combination of coherence, normalized volume and double-bounce scattering have a high correlation with SD (R² = 0.84, RMSE = 0.18). Additionally, temporal and spatial variability in SD was also observed from three polarimetric SAR images acquired between 4 April 2015 and 15 May 2015 over the Western Nordenskiöld Land region. Increase in snow depth corresponding to snow precipitation events were also detected using the POLSAR data.

1. Introduction

Seasonal snow accumulation on glaciers plays an important role as the input for mass and energy balance, as well as in controlling the thermal state of glaciers. The information on snowpack parameters is also decisive for forecasting of snowmelt runoff [1] and avalanche activity [2] in snow-covered mountainous areas. Furthermore, local meteorological [3] and topographical [4] conditions control distribution of snow accumulation and spatial variability in snowpack parameters, for example depth and density. However, field-based single-point measurements of snow depth may not be a representative of the entire distribution of snowpack accumulation over a large area. Previous investigations [5,6,7] over Svalbard have indicated that deposition of snow on the glaciers in the region is highly diverse in space due to varying topoclimatological conditions. Snow deposition increases with elevation, but it is also observed that this altitudinal trend is highly variable due to the drifting wind and other local topoclimatological factors in Svalbard. Recent investigations [8,9,10,11,12] suggest that the altitudinal gradient of snow accumulation has decreased, and the average snow thickness on the Austre Grønfjordbreen has increased by 17 cm compared to average snow thickness in 1979. Spatial variability has also increased over the region. Hence, continuous monitoring of the spatial and temporal variability of snow cover depth is the need of the hour, which can be used to understand the state of the glacier and the changes and processes taking place in it.

An alternative method to field-based manual measurements of snow depth is a ground penetrating radar (GPR) survey. This survey is an optimal technique to provide a spatial profile over large areas for understanding snow depth variability [8,9]. Apart from the spatial variability, snow accumulation and its parameters also change with time which require frequent large area observations to monitor their spatial and temporal variabilities. The importance of studies of snow cover on glaciers is not limited to only water resource assessment as the response of snow cover as an indicator of climate change is also important. Similarity/variability in the pattern and behaviour of snow cover distribution during an accumulation season(s) might be also considered in understanding of thermal state of glaciers and solid precipitation redistribution [9]. A frequent GPR survey (ground or air platform based) will not only be costly but also limited in spatial and temporal coverage. Conversely, remote sensing techniques using satellites is a cost-effective monitoring system that can be used frequently to retrieve snowpack accumulation (depth) variability in time and space [13,14,15,16,17,18,19].

Passive microwave radiometers are used to generate global snow water equivalent (SWE) and snow depth (SD) at coarse resolutions [20,21,22], but fails to capture spatial variability at the catchment and glacier-system levels. Retrieving the SD is one of the critical problems faced by the global scientific community, and many efforts have been taken in the past to develop several SD retrieval algorithms using high-resolution spaceborne active microwave synthetic aperture radar (SAR) sensor measurements. The X- and C-band SAR measurements are sensitive to snow characteristics and have been widely used to monitor snowpack density and wetness [23,24]. It also has been demonstrated that the SD can be related to the co-polar phase difference (CPD) and anisotropy factors [25,26]. But the same approach is difficult for SD inversion due to difficulties in the estimation of anisotropy. The SWE of shallow snow cover can be retrieved using Sentinel-1 SAR data [26], using the change detection technique and regression equation based on the thermal conductivity of the snowpack. However, this approach cannot be used to estimate the SD directly. The extinction coefficient calculated from the imaginary part of the dielectric constant enables the SD estimation from the X-band full polarization SAR data. However, selection of the correct model to relate the dielectric constant to wetness to retrieve the imaginary part of dielectric constant is difficult.

The retrieval of snow wetness and density based on physical approaches using fully polarimetric SAR (POLSAR) data are well established [13,14,15,16,17,18,19,23,24]. However, there is no operational SD estimation algorithm available based on POLSAR data. Some efforts have been taken in the past to relate SD with fully polarimetric SAR decomposition parameters (entropy and alpha) and the co-polar phase difference [25]. However, the relationships are not unique. Recently, a global SD mapping algorithm based on Sentinel-1 dual polarization data has been proposed [16] at 1 km spatial resolution. This method depends on cross polarization variability with snow depth, which can also be affected by the polarization orientation shift due to azimuthal slopes causing large depolarizations (cross polarization term) in mountainous regions [27]. This polarization orientation shift in snow-covered mountainous terrain may be translated into snow depth inversion. However, polarization orientation shifts can be minimized by utilizing POLSAR data and methods directly [28]. Furthermore, the POLSAR-based algorithm takes into account the local incident angle effect in the snowpack parameter retrieval process, and hence considers the local topography variation as required in the mountainous regions [13]. However, POLSAR data decomposition methods are independent of the frequency and incident angle, and they are then more suitable for snow-cover studies in mountainous terrains [14,17]. In this paper, we are proposing POLSAR-based snow depth estimation by utilizing co-polarization coherence and physical model-based scattering power decomposition.

For this purpose, we analysed L-band fully polarimetric SAR data sets acquired from the Advanced Land Observation Satellite-2 (ALOS-2) Phased Array L-band Synthetic Aperture Radar-2 (PALSAR-2) sensor over the Svalbard region. Near real-time GPR measurements over Austre Grønfjordbreen are utilized to derive relationship between the SD and the polarimetric parameters. The data is inverted to develop an empirical SD retrieval model. Further, the developed model is implemented to predict the SD and it spatial and temporal variability on other glaciers in Svalbard.

2. Study Area

As a test site we chose Austre Grønfjordbreen, located in the western part of Nordenskiöld Land (West Spitsbergen), in the head of the Grønfjorden bay (Figure 1). The length of the glacier is about 6.2 km with an average width of 1–1.2 km, with an area of 7.34 km². It ranges in elevation from 40 to 490 m a.s.l. The glacier consists of two branches merging in the middle part at 300–360 m a.s.l. The western branch of 5 km length has in its upper reaches a common ice divide with the Fridtjovbreen at 420 m a.s.l. The Austre Grønfjordbreen has been shrinking in recent decades. Between 1990 and 2011, its surface became lower by 36 m and its tongue retreated by 700 m, and by an additional 400 m between 2011 and 2015. Particularly, notable changes occurred in the region during last 10–12 years, when in summertime the snow line sometimes rose above 500 m, and the glacier changed completely in the ablation area. Now the firn zone is observed only as small patches at the ice divide and in a small area in the head of the eastern branch.

Results of previous snow surveys on the glacier (1979, 2011–2014) show that the mean snow thickness there varies in range from 140 to 165 cm (0.55 to 0.70 m water equivalent (w.e.)). Ground-based radar (GPR) measurements at 500 MHz of snow cover depth in spring of 2014 showed a close relationship (R² = 0.98) with data from manual measurements [8]. Now the firn zone is observed only as small patches at the ice divide and in a small area in the head of the eastern branch. The limited firn area was also retained in 1998 and 2000 [29]. In comparison with the data of manual snow surveys, radar measurements show a similar but more detailed picture of the distribution of snow cover depth. The discrepancy of manual and RES snow surveys in 2014 on most part of the glacier (95.7% of total area) is less than 30 cm. In April 2014 the average snow thickness on the glacier was 150.6 cm from manual surveys (77 measurement points) and 155.8 cm from the RES survey (34,754 measurement points). According to these measurements, the altitudinal gradient of snow accumulation at Austre Grønfjordbreen in 2014 was 0.21 m/100 m, which is less than its mean value there (0.35 m/100 m). In addition to an altitudinal gradient, there is asymmetry in snow distribution relative to the long axis of the glacier. Generally, more snow accumulates in the western part than in the eastern one, which seems to reflect the pattern of local wind regime, snow drift, and deposition.

3. Data Used

3.1. GPR Data and Snow Measurements

In this study, we use snow measurements carried out between 7 and 13 April 2015, before snow melting started, when the entire thickness of the snow cover had negative temperatures. The standard snow survey measurement points consisted of measurements of snow depth using a snow probe and were carried out along the axis of the glacier in its central part, and their coordinates were determined using a GPS receiver with an accuracy of 5 m. Additionally, thickness and average density of the snow cover were measured in 12 pits located in the lower, middle, and upper part of the glacier, with measuring the mass of snow samples taken layer-by-layer using a tube sampler of 100 cm³ volume.

For the radar snow thickness survey performed on 10 April 2015, we used the pulseEKKO PRO GPR with 500 MHz antennas. GPR includes a receiving and transmitting device with power supplies of 12 volts, receiving and transmitting shielded antennas, a control unit with digital registration of radar data, and a GPS receiver mounted on a wooden sledge transported by a snowmobile. Starting the transmitting device is carried out either by a timer with an interval of 1–2 s or from the odometer with an interval of about 1.3 m at a speed of 5–10 km/h on the glacier. At a frequency of 500 MHz (wavelength in the air λ = 0.6 m), the transmitter generates electromagnetic broadband pulses of 400 volts, duration 1 ns with a repetition rate of 100 kHz, which are reflected from the dielectric contrast interfaces in the snow cover, received by the receiver in a time window of 50 ns and digitized with a sampling period of 0.2 ns. At an average radio wave propagation speed of 0.2 m/ns, this allows to obtain reflections from the boundaries in the snow cover to a depth of 6.25 m with a maximum vertical resolution λ/4 = 0.1 m and a limiting resolution λ/0 = 0.04 m. Measurements along the profiles were performed at a fixed distance between the centres of the transmitting and receiving antennas of x₀ = 0.23 m. For visualization and interpretation of radar reflections, we used the software package RadexPro Plus 2011.2 Basic 2011.1 [30]. The total length of the snow survey profiles on the glacier was near 44 km with 36,794 points of measurements (see Figure 1).

To convert two-way travel times of reflected radar waves to snow depths we used our data of snow density measurements in 12 snow pits and Looyenga’s formula [31] for heterogeneous mixtures that gives the best estimation compared with other solutions as we found earlier [8]. The averaged mean snow density in snow pits varies from 359 to 580 kg m⁻³. Highest density was registered in two pits (located at the ice divide and terminus) where snow cover contains 3–5 ice layers (with mean thickness 2.5 cm each). For time-to-depth conversion we chose the mean value of density from all 12 pits (433 kg m⁻³) that corresponds to the radio wave velocity in snow at 0.221 m/ns.

3.2. ALOS-2/PALSAR-2 Data

The state-of-the-art PALSAR-2 aboard ALOS-2 is an active microwave radar operating at 1.2 GHz of the L-band frequency range, which have enhanced performance compared to its predecessor ALOS/PALSAR and JERS-1/SAR [32]. Launched on 24 May 2014, ALOS-2/PALSAR-2 operates at an altitude of 628 km with a revisit time of 14 days. Expanding its transmission power and bandwidth and along with adopting newer techniques such as dual-beam receivers, complex chirp modulations, and highly-efficient data compression, ALOS-2/PALSAR-2 provides data with higher resolution, wider swath, and improved image quality.

As per the requirements of the earth observation application, PALSAR-2 operates in different observation modes of spotlight, stripmap, scanSAR, and full polarimetry mode. The observation modes vary in polarization, swath width and spatial resolution in range, azimuth directions, and bandwidth. In our study, we have utilized fully polarimetric data acquired in the stripmap mode at a resolution of 5.1 m in range direction and 4.3 m in azimuth direction. The details of acquired data utilized in this study are listed below in

Table 1. The details of acquired ALOS-2/PALSAR-2 data sets.

Date of Acquisition	Polarization	Off-Nadir Angle	Range Resolution	Azimuth Resolution	Level of Data Product	Observation Mode/Orbit Pass
4 April 2015	Quad	32.7°	5.1 m	4.3 m	SLC, Level 1.1	Stripmap/Ascending
13 April 2015	Quad	30.4°
15 May 2015	Quad	25.0°

.

4. Methods

4.1. ALOS-2/PALSAR-2 Data Processing

In this paper, fully polarimetric ALOS-2/PASAR-2 data, acquired on 4 April, 13 April, and 15 May 2015, were used to analyse the snow depth (SD) behaviour with polarimetric parameters over the Svalbard region. Available as single-look–complex (SLC), the level 1.1 data was radiometrically calibrated. The calibration of SLC data was performed for all four elements of the 2 × 2 scattering matrix.

The 2 × 2 scattering matrix [S] obtained as the initial radar measurement from the data was not suitable for naturally occurring distributed targets, like snow and glaciers [33,34]. Hence, [S] was converted into a second-order relative phase coherency matrix [T] using Pauli ($k$) target vectors. Elements of these matrices are further used to define the polarimetric descriptors/parameters with which relationships with SD are established.

• Data calibration and matrices formation

The complete target scattering information is defined by the relation between the incident ($\overrightarrow{E_I}$) and scattered ($\overrightarrow{E_S}$) electromagnetic (EM) wave fields, which is given by the scattering matrix [S] for HV linear polarization basis as [35]

$$E_S=\frac{e^{-jkr}}{r}\left[S\right]E_I=\frac{e^{-jkr}}{r}\left[\begin{array}{cc}{S}_{HH}& S_{HV}\\ S_{VH}& S_{VV}\end{array}\right]E_I.$$

(1)

Available as single-look–complex (SLC) data from ALOS-2/PALSAR-2, the data is first calibrated. The calibration of the SLC [S] data matrix can be generated using the following equation:

$${\left[S\right]}_{cal}=k\left[\begin{array}{cc}{S}_{HH}& S_{HV}\\ S_{VH}& S_{VV}\end{array}\right] ,$$

(2)

where, $k=\sqrt{{10}^{\left(CF-32\right)/10}}$ is the factor used for radiometric calibration of the data. The calibration factor CF can be found from the header file of the data.

Furthermore, the second-ordered scattering matrix [35] can be formed from ${\left[S\right]}_{cal}$ based on Pauli basis matrices by considering spatial ensemble averaging of a factor of 4 in range and 6 in azimuth directions respectively as

$$\begin{array}{c}\langle \left[T\right]\rangle =\langle k{k}^{†}\rangle =\left[ \begin{array}{ccc}\langle {\left|k_1\right|}^2\rangle & \langle k_1{k}_2{}^{*}\rangle & \langle k_1{k}_3{}^{*}\rangle \\ \langle k_2{k}_1{}^{*}\rangle & \langle {\left|k_2\right|}^2\rangle & \langle k_2{k}_3{}^{*}\rangle \\ \langle k_3{k}_1{}^{*}\rangle & \langle k_3{k}_2{}^{*}\rangle & \langle {\left|k_3\right|}^2\rangle \end{array}\right]=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right]=\\ \frac{1}{2}\left[\begin{array}{ccc}\langle |S_{HH}+S_{VV}{|}^2\rangle & \langle (S_{HH}+S_{VV}){(S_{HH}-S_{VV})}^{*}\rangle & 2\langle (S_{HH}+S_{VV})S_{HV}{}^{*}\rangle \\ \langle (S_{HH}-S_{VV}){(S_{HH}+S_{VV})}^{*}\rangle & \langle |S_{HH}-S_{VV}{|}^2\rangle & 2\langle (S_{HH}-S_{VV})S_{HV}{}^{*}\rangle \\ 2\langle S_{HV}{(S_{HH}+S_{VV})}^{*}\rangle & 2\langle S_{HV}{(S_{HH}-S_{VV})}^{*}\rangle & 4\langle |S_{HV}{|}^2\rangle \end{array}\right],\end{array}$$

(3)

where

$$k=\frac{1}{\sqrt{2}}\left[\begin{array}{c}{S}_{HH}+S_{VV}\\ S_{HH}-S_{VV}\\ 2S_{HV}\end{array}\right]=\left[\begin{array}{c}{k}_1\\ k_2\\ k_3\end{array}\right] ,$$

(4)

and in which * indicates the complex conjugation, † denotes the complex conjugation and transposition, and 〈…〉 indicates the spatial ensemble averaging which is represented here by multilook factors in the azimuth and range directions.

b.

Orientation angle compensation in [T]

Compensation of the rotation angle applied on [T] further reduces the number of independent [T] parameters from nine to eight. $\langle \left[T\left(\theta \right)\right]\rangle$ is the [T] matrix rotated about the radar line of sight, and is given as [28]

$$\langle \left[T\left(\theta \right)\right]\rangle =\left[\begin{array}{ccc}{T}_{11}\left(\theta \right)& T_{12}\left(\theta \right)& T_{13}\left(\theta \right)\\ T_{21}\left(\theta \right)& T_{22}\left(\theta \right)& T_{23}\left(\theta \right)\\ T_{31}\left(\theta \right)& T_{32}\left(\theta \right)& T_{33}\left(\theta \right)\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos 2\theta & \sin 2\theta \\ 0& -\sin 2\theta & \cos 2\theta \end{array}\right]\langle \left[T\right]\rangle \left[\begin{array}{ccc}1& 0& 0\\ 0& \cos 2\theta & -\sin 2\theta \\ 0& \sin 2\theta & \cos 2\theta \end{array}\right] ,$$

(5)

where

$$2\theta =\frac{1}{2}{\tan }^{-1}\left(\frac{2\Re \left\{T_{23}\right\}}{T_{22}-T_{33}}\right).$$

(6)

This is followed by implementation of Refined Lee filter, with a window size of 5 × 5 pixels, to further reduce the speckle noise and then a terrain correction is performed on the coherency matrix by utilizing a 30 m posting ASTER-Global Digital Elevation Model (ASTER-GDEM).

c.

Six-component scattering matrix power decomposition (6SD)

The POLSAR decomposition methods in radar polarimetry are classified into two categories: (1) coherent and (2) incoherent [33]. The coherent target decomposition expresses the measured scattering matrix in terms of contribution from simple canonical objects (coherent objects). However, it is used to characterize the [S] matrix only. The incoherent decomposition uses a second-order polarimetric matrix-like coherency matrix [T] rather than an [S] matrix, and it is applicable to distributed targets (natural objects) [34]. The incoherent model-based SAR decomposition theory is well established and started with Freeman Durden Decomposition (FDD) [36], which decomposes the total power received by the SAR sensor into simple canonical targets (surface, double-bounce, and volume scattering). The surface scattering (single/odd-bounce scattering) is caused by smooth surfaces or first-order Bragg scatterers. The dihedral scattering (even/double-bounce scattering) is observed from corner reflectors [37], e.g., between the road and wall of the building. In double-bounce scattering, the incident electromagnetic (EM) wave encounters two reflections resulting in a 180° shift between the incident and reflected EM wave. When an EM wave interacts with complex structures like trees and snow, the multiple reflections within the medium, due to dielectric discontinuity, lead to the volume scattering.

In case of a snow or dry snow-covered glacier, the surface scattering can be observed from air-snow, or the snow-ice interface [14]. On the other hand, the volume scattering in snowpack occurs due to multiple reflections of EM waves. However, double-bounce scattering is small as compared to the surface and volume scattering in case of snow [13,14], but in a glaciated region, the crevasse and extremely rough interfaces may lead to significant double-bounce scattering contribution, because microwave radar signals passes through dry snowpack and interacts with the glacier surface [27,37].

Utilization of model-based decompositions for this study is due to their ease of implementation followed by an accurate interpretation and classification of naturally occurring targets and physical relevance of scattering powers to the scatterers. Furthermore, the advancement in model-based scattering power decompositions [28,38,39,40,41] are achieved, and for the purpose of this study, we have adopted the latest six-component scattering matrix power decomposition of the coherency matrix $\langle \left[T\left(\theta \right)\right]\rangle$ to generate six scattering powers from POLSAR data. The six-component scattering power decomposition model (6SD) separates six scattering powers from the composite total power (TP) [41] as given below:

$$\begin{array}{c}\langle \left[T\left(\theta \right)\right]\rangle =P_{S }{\left[T\right]}_S+P_D{\left[T\right]}_D+\{\begin{array}{c}{P}_V{\left[T\right]}_{vol}^{dipole} \\ P_V{\left[T\right]}_{vol}^{dihedral} \end{array}+P_H{\left[T\right]}_H+P_{OD}{\left[T\right]}_{OD}+P_{CD}{\left[T\right]}_{CD}\\ =\frac{P_{S }}{1+|\beta {|}^2}\left[\begin{array}{ccc}1& {\beta }^{*}& 0\\ \beta & |\beta {|}^2& 0\\ 0& 0& 0\end{array}\right]+\frac{P_D}{1+|\alpha {|}^2}\left[\begin{array}{ccc}|\alpha {|}^2& \alpha & 0\\ {\alpha }^{*}& 1& 0\\ 0& 0& 0\end{array}\right]\\ +P_V\left\{\begin{array}{c}\frac{1}{4}\left[\begin{array}{ccc}2& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]; \mathrm{random} \mathrm{distribution}\\ \frac{1}{30}\left[\begin{array}{ccc}15& -5& 0\\ -5& 7& 0\\ 0& 0& 8\end{array}\right]; \mathrm{horizontal} \mathrm{distribution}\\ \frac{1}{30}\left[\begin{array}{ccc}15& 5& 0\\ 5& 7& 0\\ 0& 0& 8\end{array}\right]; \mathrm{vertical} \mathrm{distribution}\\ \frac{1}{15}\left[\begin{array}{ccc}0& 0& 0\\ 0& 7& 0\\ 0& 0& 8\end{array}\right]; \mathrm{oriented} \mathrm{dihedral}\end{array}\right\}\\ +P_H\frac{1}{2}\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& \pm j\\ 0& \mp j& -1\end{array}\right]+P_{OD}\frac{1}{2}\left[\begin{array}{ccc}1& 0& \pm 1\\ 0& 0& 0\\ \pm 1& 0& 1\end{array}\right]+P_{CD}\frac{1}{2}\left[\begin{array}{ccc}1& 0& \pm j\\ 0& 0& 0\\ \mp j& 0& 1\end{array}\right] ,\end{array}$$

(7)

where the six scattering powers, P_S, P_D, P_V, P_H, P_CD, and P_OD, correspond to surface scattering, double-bounce scattering, volume scattering, helix scattering, compound dipole scattering (oriented quarter-wave reflector scattering), and oriented dipole scattering, respectively. The solution to the above equations and details about the expansion of the individual scattering submatrices, ${\left[T\right]}_S$, ${\left[T\right]}_D$, ${\left[T\right]}_{vol}$, ${\left[T\right]}_H$, ${\left[T\right]}_{OD}$, and ${\left[T\right]}_{CD}$, are given in [41].

Using the six powers (surface (P_S), volume (P_V), double-bounce (P_D), helix (P_H), oriented dipole (P_OD), and compound dipole (P_CD)), three polarimetric parameters were derived, namely the normalized volume scattering power (P_NV = P_V/TP), normalized double-bounce scattering power (P_ND = P_D/TP), and volume to double-bounce scattering ratio (P_VD), which were used to see the potential of model-based decomposition for SD retrieval.

d.

Coherence

The polarimetric coherence ($\overline{{\gamma }_c}$) was also calculated using co-polarization information as in [35]:

$$\overline{{\gamma }_c}={\gamma }_c{e}^{i{\phi }_c}=\frac{S_{HH}{S}_{VV}^{*}}{\sqrt{S_{HH}{S}_{HH}^{*}}\sqrt{S_{VV}{S}_{VV}^{*}}} ,$$

(8)

where $\overline{{\gamma }_c}$ is complex coherence, ${\gamma }_c$ is absolute coherence, and ${\phi }_{c }$ is co-polar phase difference. ${\gamma }_c$ (denoted by COH (${\gamma }_{HH,VV}$) hereafter) is utilized for further analysis to establish the SD relation and retrieval.

4.2. Field-Measured GPR Data Processing and Data Points Reduction

An extensive field campaign was conducted over the test site in April 2015 as mentioned in Section 3.1. The SD was measured at a total of 36,794 points taken at a sampling interval of less than one meter. Due to the multilooking and spatial resolution of the SAR data, one SAR pixel covers several GPR measurements on the ground. In an effort to bring the measurements from SAR and GPR onto a common spatial scale, an averaging of GPR points is done corresponding to the spatial size of a single SAR image pixel (30 m). After applying the averaging procedure on GPR measurements corresponding to the same SAR resolution cell, the total ground measured points were reduced to 2288, denoted as G (Figure 2a).

Further, the G data set is divided into two subsets using alternate sampling for testing and validation [42]. These subsets are denoted as G1 and G2, each with 1144 samples (N), shown in Figure 2b,c, respectively. The histograms of G1 and G2 are shown in Figure 2d,e, respectively. G1 and G2 have a similar SD distribution with approximately the same mean, minimum, and maximum value ranges, indicating unbiased sampling of the GPR points. To capture the heterogeneity of the measured SD with respect to polarimetric SAR descriptors, G1 and G2 were averaged by dividing polarimetric coherence (range 0 to 1) into 100 classes with 0.01 increments in coherence. For example, class 1 varies between 0 to 0.01, class 2 covers a range of coherence values from 0.01 to 0.02, and so on. This reduction strategy in GPR points is important and adopted for developing relations between SD and polarimetric parameters as discussed in Section 5.3.

5. Results

5.1. 6SD Interpretation

A six-component scattering matrix power decomposition (6SD) was applied on ALOS-2/PALSAR-2 data sets. Figure 3 shows the RGB image, formed by assigning a standard colour code to the model-based decomposition powers (red: double-bounce scattering; green: volume scattering; and blue: surface scattering). The colour-coded image indicates mostly volume scattering dominance over the glaciers in Figure 3. Dominant double-bounce scattering is observed over a few points, such as the crevasses zone of the glaciers for example the frontal areas of Vestre Grønfjordbreen, Fridtjovbreen, and other glaciers. L-band radar signal penetrating through the snowpack can significantly interact with crevasses, which produces dihedral type scattering. The blue colour pattern is also observed on Austre Grønfjordbreen near the terminus, which detects dominate surface scattering contribution in this part of glacier. A special pattern of cyan and white colours is noticed clearly in the firn area at the head (higher elevation) of the glaciers. It is formed due to strong contributions of surface, volume, and double-bounce scattering powers in the firn region of glacier. In general, dry snow over ground is generally transparent and is not deep enough to give high returns at the L-band. But a thick dry snow and firn in the accumulation areas of glaciers provide higher backscattering because of low absorption of the signal and the likelihood of scattering by large snow grains (the Rayleigh scattering regime of microwaves). It is observed that Vestre Dahlfonna, Erdmannbreen, Fridtjovbreen, and Vestre Grønfjordbreen exhibit a noticeable portion of firn area but the firn area ratio is small as compared to the total area of these glaciers. Firn area ratios of Austre Grønfjordbreen and Austre Dahlfonna are also not considerable. However, the variation in the surface, volume, and double-bounce scattering noticed over Austre Grønfjordbreen is similar to other glaciers, as shown in Figure 3. To get the understanding and the spatial variation of all six scattering powers, the map of individual scattering powers (surface, volume, double-bounce, oriented dipole, complex dipole, and helix) over Austre Grønfjordbreen are shown in Figure 4. Variation in compound scattering powers (helix, oriented dipole, and compound dipoles) represents heterogeneity of the medium. It can also be related to spatial variation in deposition of snow over the glacier.

Furthermore, the variations of scattering powers are studied in detail to understand their behaviours from the terminus toward head (see Figure 4) on Austre Grønfjordbreen. In Figure 4, P_H, P_OD, and P_CD show small-scale random variations. These are the minor components of the 6SD model decomposition, which dictates the heterogeneity in the medium. If these components are zero, the medium is said to be homogeneous. The behaviour of 6SD decomposed powers along the transact from A to B is also studied with the help of power profiles. Figure 5 shows the variation of normalized powers along transact (A to B, see Figure 3). At the beginning of the transact (near the terminus), the surface scattering is dominating over volume and double bounce. As the snow depth starts increasing along the transact, the volume scattering starts to dominate over other scattering powers. There is not much variation in the double-bounce scattering but it shows an overall decreasing trend along the transect from position A to B. These variations of scattering powers in snowpack over the glacier can be associated with variations in accumulation rate and layer density and grain size (growth with age) [43,44].

The firn is characterized by a vertical structure and hence shows a large anisotropy effect. A very negligible change is observed in helix scattering (P_H). The trend line also suggests that oriented (P_OD) and compound dipole scattering (P_CD) increases towards point B of the transect. Oriented dipole (P_OD) and compound dipole scattering (P_CD) also shows large fluctuations in the firn region (towards point B). These fluctuations can be related to an inhomogeneous accumulation rate and associated anisotropy effect. The inhomogeneous accumulation and settling lead to the horizontal and vertical orientation of the snow particles in the snowpack [25]. This leads to dielectric anisotropy and fluctuating returns observed as oriented and compound-oriented dipole scattering.

5.2. Coherence Image

Polarimetric coherence over the western part of Nordenskiöld Land is shown in Figure 6. The image shows the spatial variability of coherence along the glaciers. Six glacier boundaries are shown with Austre Grønfjordbreen (marked by a black outline) shown separately in Figure 6 (right).

Polarimetric coherence represents the similarity of the scattering measured at the two polarimetric channels (HH and VV). Targets/features with canonical structures have lower polarimetric coherence as the two polarimetric channels will observe distinct scattering pattern. Over snow, as the snow depth increases, the complexity of the scattering mechanism also increases, leading to similar scattering from HH and VV polarized channels. This leads to increased polarimetric coherence. Thus, there is a direct correlation between polarimetric coherence and snow depth. However, an exception is to be made for specular scatterers. For such surfaces, there is very low backscattering towards the sensor for both the HH and VV polarizations. Although the backscattering is low, the polarimetric coherence can be high. For the Austre Grønfjordbreen, the coherence is higher at the terminus (Figure 6). Furthermore, as observed in Section 5.1, this terminus region shows high surface scattering. The mean backscattering, ${\gamma }_{HH}$ and ${\gamma }_{VV}$, for the terminus are −14.17 dB and −13.59 dB, respectively. At the terminus, the ground measured snow depth is low (<0.75 m), exposing the incident microwaves to the ice surface beneath the snow layer. This leads to specular scattering, high surface scattering, low backscattering, and high polarimetric coherence. For such regions, the polarimetric coherence is not correlated with snow depth.

5.3. Behaviour of Polarimetric Parameters with Snow Depth

We have identified four polarimetric parameters, namely polarimetric coherence COH (${\gamma }_{HH,VV}$), ratio of the volume to double-bounce scattering powers ($P_{VD}$), normalized double-bounce scattering power ($P_{ND}$), and normalized volume scattering power ($P_{NV}$), to study their relationship with snow depth. For this purpose, G1 and G2 sets of GPR measurements have been utilized. The relationships between polarimetric parameters and SD for the G1 set of points are shown in Figure 7.

A positive correlation between SD with COH and P_NV can be seen in Figure 7a,d, respectively. With the increase in SD, the P_VN component increases due to the increased volume scattering from the snowpack. It is observed that the logarithmic fit between SD and $P_{VD}$ (shown in Figure 7b) gives good correlation, explaining the SD behaviour with both volume and double-bounce scattering power. At the same time, the P_ND component decreases with SD due to less contribution of double-bounce scattering, as shown in Figure 7c. As the snow depth increases, the volume scattering contribution to the total scattering increases, at the same time the surface and double-bounce scattering decreases as shown Figure 5. Volume scattering in the snow increases with the degree of heterogeneity in snowpack and its depth. This thick heterogeneous snowpack (i.e., multi-layered snow deposit varying in terms of its bulk properties) can increase the multiple bounce scattering from the snowpack and lead to depolarization phenomena due random processes of scattering. These effects can decrease single- and double-bounce scattering with increasing heterogeneous snowpack depth. The surface and double-bounce scattering show a negative correlation with SD across the test glacier, for example Figure 7c depicts a double-bounce trend with SD. As COH demonstrates a strong positive correlation with SD, it indicates that the difference between HH and VV backscattering reduces due to similar backscattering at higher SD. Hence, the polarimetric coherence between HH and VV polarization increases with SD. As explained in Section 5.2, the COH is a strong indicator of SD. The behaviour of SD with polarimetric parameters was also tested for the G2 set of GPR data and it is shown in Figure 8. Similar trends are observed for the G2 set of GPR data.

5.4. Inversion and Validation

The earlier section analysed the relationship between various polarimetric parameters and SD. Using this analysis, SD inversion algorithms are developed, trained, and validated in this section. Figure 7 and Figure 8 show four polarimetric parameters—COH-, P_V/D-, P_NV-, and P_ND-based relationships with SD. For the development of the SD inversion algorithm, the relationships shown in Figure 8a is adopted initially. All other possible combinations of polarimetric parameters were evaluated and a few selected models, which showed consistent performance, are discussed here. While the validation of SD inversion using individual-parameter-based relations is tested, we utilize a two-parameter-based model with inputs of normalized double-bounce $P_{ND}$ and normalized volume scattering $P_{NV}$. Further, one model with three polarimetric parameter inputs of COH, $P_{ND}$, and $P_{NV}$ is also analysed. In total, we analysed six regression-based models for data sets G1 and G2 independently. As explained in Section 4.2, the GPR measured SD and polarimetric descriptors are divided into two sub-sets of data (i.e., G1 and G2). These data sets are used as training and validation sets for analysis of the six models. The regression models are trained using the G1 data set and validated with the G2 data set. Additionally, to test the stability of the developed models, the training and validation data are reversed, and the models are trained using GPR data in G2 and validated on data in G1. The model inverted SD is validated with GPR-measured SD by two parameters: root mean squared error (RMSE) and coefficient of determination (R²). The results for the two training and validation cases (G1 and G2) for all the six models is shown in

Table 2. Validation results of regression analysis-based SD inversion models. The R² and RMSE of the six models is given for two cases: training with the G1 set and with the G2 set.

Parameters		COH	${\mathit{P}}_{\mathit{N}\mathit{D}}$	${\mathit{P}}_{\mathit{N}\mathit{V}}$	${\mathit{P}}_{\mathit{V}/\mathit{D}}$	${\mathit{P}}_{\mathit{N}\mathit{D}},{\mathit{P}}_{\mathit{N}\mathit{V}}$	$\mathit{C}\mathit{O}\mathit{H},{\mathit{P}}_{\mathit{N}\mathit{D},}{\mathit{P}}_{\mathit{N}\mathit{V}}$
Trained G1 Validated G2	R2	0.84	0.72	0.73	0.83	0.84	0.84
	RMSE (m)	0.18	0.23	0.23	0.18	0.18	0.18
Trained G2 Validate G1	R2	0.78	0.78	0.64	0.78	0.82	0.82
	RMSE (m)	0.21	0.20	0.26	0.21	0.19	0.19

.

Figure 9 show the validation plots of GPR-measured SD with model inverted SD for all six models. In Figure 9, the GPR data in G1 is used to train the models, while the model inverted SD is validated using the G2 data set (shown in black colour). Similarly, the GPR data in G2 is used to train the models and the data set in G1 is used for validation (shown in red colour). The polarimetric coherence, ratio of volume to double-bounce scattering powers, and normalized double-bounce power demonstrate to be reasonably good parameters to estimate SD in Figure 9. However, the inversion accuracy is found to be low with the normalized volume scattering power, because variation in the normalized volume scattering power between 0.3 to 0.4 is less sensitive with a shallow depth snowpack.

To demonstrate multivariate SD inversion, we have used the combinations of polarimetric parameters as shown in Figure 9. SD inversion results from the combinations of P_ND and P_NV shown in Figure 9e provides a better estimate of SD compared to the individual polarimetric parameters. The other two parameter combination models still give good accuracy. But the combination of normalized volume and normalized double-bounce scattering power provided the best accuracy with low RMSE. The range of the P_ND parameter is restricted to 0.4, and P_NV shows variation between 0 and 1.

Hence, inversion based on the combination of P_ND and P_NV do not cover the full range (0 to 1) of P_ND variations, leading to negative SD values at some available ground measurements locations. To avoid this problem, we have added a third variable of COH to an existing two-variable model. This combination of polarimetric parameters leads to the non-negative SD inversion over the entire study region but the inclusion of the third parameter does not improve the SD accuracy. Similar behaviours of developed regression model based on another set of data (i.e., regression equation of the G1 data set is used to invert SD from G2 data set) are observed. Hence, polarimetric scattering power-dependent relations are not used for further analysis in study. Apparently, scattering decomposition-dependent relations have shown potential to develop an SD inversion model (which need to be investigated further in a separate study) but the magnitude of the COH-based regression model illustrates a well-linearized function of varying SD. Therefore, further analysis was achieved based on the COH magnitude and SD relation.

5.5. Comparison SAR and GPR-Based SD and Differences

After analysing all the combination for the SD estimation. We have chosen one of the simple and best regression relations, which is based on co-polarization coherence. It is given as

$$\mathrm{SD}=\left(2.2006\times COH\right)+0.5661.$$

(9)

Generally, the SD inversion results from POLSAR and interpolated RES snow survey data are in fairly good agreement (Figure 10). Of the total glacier area, 36% area has less than 30 cm discrepancy between RES (13 April 2015) and POLSAR snow thickness on 10 April 2015. Likewise, 59% of the area has discrepancy less than 50 cm. Mean snow thickness from RES interpolation is 1.49 m (range between 0.23 and 3.68 m), and from SD inversion is 1.8 m (range between 0.57 and 2.74 m).

SD inversion shows higher values than measured. In addition, difference maps were generated (Figure 11) to compare results of SD from POLSAR data with interpolated SD from GPR. Figure 11 demonstrates a large positive difference at the terminus where shallow depth snow accumulation exists. It is also observed that volume and double-bounce scattering are low and surface scattering is higher in the same region. Low volume scattering represents that there are no random scattering processes that increase with SD in general. For the SD of 1.75 m or above, the radar backscatter of HH and VV polarizations is dominantly coming from the snowpack (volume scattering contribution >50%) with a similar strength of amplitude, resulting in the high coherence value (>0.5). For the SD range of 0.75 to 1.75 m, the HH and VV backscattering are contributed from different layers within the snowpack and ice. Volume scattering contribution varies between 20% and 50%. Similar variations, also noted in POLSAR coherence, have a linear relationship with SD ranging between 0.57 m and 2.74 m (Figure 7).

However, for the snow thickness <0.75 m, L-band microwave SAR signals penetrate through the snowpack (volume scattering contribution is negligible) and interact at the snow–ice interface, leading to scattering similar to specular reflections and low radar signal return amplitudes, which is similar in both HH and VV polarizations (high coherence). High POLSAR coherence quantifies the similarity between the backscattering measurements at HH and VV polarizations. Near the terminus snow depth is low, and the L-band radar backscattering contribution may come from the same/similar physical structure beneath thin snow layer. This cause is producing high POLSAR coherence, which is translating into SD inversion results on the higher side. The proposed coherence-based algorithm performs poorly in areas with shallow snow depth over the glacier (such as the terminus where the scattering contribution from the snowpack is less), leading to a high discrepancy in SD inversion results.

Furthermore, the maps showing the difference in snow thickness measured from POLSAR data and the RES points, along with the error maps for each SAR acquisition, also indicate negative errors in the south-west part and at the southernmost tip. A negative difference (i.e., underestimation) in these regions might be due to the uncertainty of field data extrapolation in the area not covered by RES. However, the difference is the smallest in this region along the GPR data profile. Therefore, SD inversion for this area might be closer to reality than the extrapolated field data. With some uncertainty, the proposed methodology based on SAR data gives an alternate approach for SD estimation, which is more feasible than manual and GPR measurements. The spots demonstrate that errors of underestimation occur near those stretches where the sounding profiles observed large spatial differences in snow thickness cover compared to its surroundings. The resulting map (see Figure 10) dictates that it is important to interpret how optimally inverted SD from POLSAR data varies along the GPR profile.

5.6. Spatial and Temporal Variability in Snowpack Depth

The SD inversion maps over the glaciers in the western part of the Nordenskiöld Land region and the test region using Equation (9) is shown in Figure 12. In addition to the snow depth estimated, we also demonstrate the spatial and temporal variability in SD. This is performed using POLSAR data acquired over the three dates of 4 April, 13 April, and 15 May in 2015. Figure 12a–c show the inverted SD maps for the western part of Nordenskiöld Land region, Svalbard area on the three dates. Between 13 April 2015 and 15 May 2015, several significant snow events were recorded. The total amount of snow precipitated between 2 and 13 April 2015, was 10.2 mm compared to a significantly larger value of 64.6 mm of snow between 13 April 2015 and 15 May 2015. Concentrating our study, on the test glacier only, we observed that the SD has increased over more regions, indicated by more snow depth (shade of blue) over the glacier. However, the area situated on the left side of the terminus shows increased snow depth in all three maps. The reason for this may be explained as follows: The presence of a big surface channel with depth of several meters, situated on the left side of terminus, assists in the accumulation of much more snow than compared to the surrounding flat glacier surface. But snow depth near the snout is overestimated due to the snow–ice interface scattering, which produces high co-polarization coherence. Furthermore, low variability in the range of SD (maximum and minimum values) for the April and May acquisitions may be explained as a limitation of the inversion model. Nevertheless, maps in Figure 12 are indicators that the proposed inversion equation is significantly capable of capturing both temporal and spatial variability of changes in SD over the Svalbard region, and especially over the area of interest of this investigation.

6. Conclusions

In this paper, the behaviour of SD with several polarimetric parameters at the L-band is studied. Fully polarimetric SAR data from ALOS-2/PALSAR-2 was used to generate the polarimetric coherence and 6SD decomposition powers. These parameters are studied with the help of GPR measurements that are carried out in near real-time. The resolution of GPR measurements are reduced to match the resolution of the SAR products. After relating the SD with the polarimetric parameters, some interesting results are observed and discussed in detail in this paper.

Polarimetric scattering powers over different areas of the glaciers help identify the different features, such as crevasse and firn areas, owing to the different dominant scattering observed for them. In particular, the firn areas demonstrate variations in scattering powers, especially in the compound scattering powers, which is indicative of the heterogeneity of the medium and the associated anisotropy effect. Furthermore, four different polarimetric parameters have been introduced to show their relationship with snow depth. Among them, the co-polarization coherence (between HH and VV) and the ratio of volume to double-bounce scattering powers show a good correlation with the snow depth. A negative correlation is found with the double-bounce scattering power. This relationship is tested with two sets of GPR points, i.e., the ones used for testing as well as the set used for validation.

The study is taken further by introducing a snow depth inversion algorithm using the combination of these polarimetric parameters. Among the univariate models for SD inversion, inputs with polarimetric coherence performed the best while the combination of coherence with normalized volume and double-bounce scattering powers demonstrated the power of the multivariate SD inversion process. An R² of 0.84 and RMSE of 0.18 are obtained with these inversions.

The inverted SD using POLSAR data are also validated with the ground-based GPR measurements, showing a fairly good agreement. The error map generated between the POLSAR inverted SD and GPR measurements indicated positive errors near the terminus and negative errors to the south-west part of the Austre Grønfjordbreen. While shallow depth snow accumulation explains the positive error difference, uncertainty in the extrapolation of the field points could provide a plausible explanation for the negative errors. We have even extended the analysis to demonstrate the spatial and temporal variability in snowpack depth estimated over the entire Western Nordenskiöld Land region. This is performed using multi-temporal POLSAR data that span over three different dates with varied snow precipitation conditions. Between 13 April 2015 and 15 May 2015, four to five significant snow spells were recorded. Due to this, the SD on 15 May 2015 is more than on 13 April. The increase in snow depth due to some heavy snow precipitation events is well captured by the proposed inversion equation using the POLSAR data. The proposed methodology is useful in continuous monitoring of the spatial and temporal variability of snow cover depth. However, it is noticed that when the snow thickness is less than 0.75 m near the terminus, the discrepancy between POLSAR and RES interpolated snow thickness (depth) is more than 1.2 m, covering 8.7% of total glacier area. This issue with the proposed method related to shallow snow depth retrieval can be addressed in the near future.