UVicSPACE: Research & Learning Repository

Landfast


Introduction
The spatial distribution, evolution, and extent of melt ponds on sea ice are functions of snow thickness and snow redistribution processes through the preceding winter and spring, ice surface roughness, and ice type [1][2][3][4].For spatially varying snow covers on sea ice, thicker snow takes longer to melt than thinner snow under consistent atmospheric forcing [1,5,6].Prior studies advocate that a thinner snow cover leads to dominant surface flooding (i.e., larger areal melt pond fraction, f p ), whereas a thick winter snow cover leads to a greater fraction of snow patches and thus less surface flooding (i.e., smaller f p ) [7,8].Thicker snow covers generally accumulate on comparatively rough first-year sea ice (FYI), as the uneven surface topography effectively captures wind-blown snow [9].The relatively smooth topography of FYI causes thinner snow covers, with consistent drifting of snow following the wind direction during depositional storm events [10].Therefore, smooth FYI, with an associated thin snow cover, leads to the horizontal spreading of ponds over a larger area, with f p as high as 80% [3,4,11].Since the formation of FYI melt ponds is influenced by winter snow thickness, f p could provide a proxy for winter ice surface roughness and snow thickness.
After melt onset, melt ponds gradually develop over four successive stages.These stages have been previously identified as: (1) topographic control, (2) hydrostatic balance, (3) ice freeboard control and (4) fall freeze up or ice break up [6,8,12].During stage 1, melt ponds start to form and eventually fill topographic depressions with minimal drainage through the ice.At this stage, the structure of the melt ponds is very unstable and is actually controlled by the pre-melt ice topography [1].Stage 2 begins when melt ponds have the same production and drainage rates, and lateral melt water flows towards macroscopic fractures and seal holes promote the formation of melt pond networks by interconnecting melt ponds.Consequently, the geometric structure of melt pond covers becomes prominent during this stage [1,12].Stage 3 is characterized by increased vertical drainage and decreased horizontal discharge to macroscopic fractures.Surface topography of melting and decaying ice cover controls the melt pond coverage in this stage [6].Finally, during stage 4, complete decay or disintegration of FYI, or fall freeze up of multiyear sea ice begins.
Estimating melt pond fraction over large areas requires satellite-based methods.Optical approaches [13][14][15][16][17][18] are hampered by cloud cover and are thus limited by the accuracy of cloud identification and screening techniques.They must also contend with regional differences in spectral properties [18], and issues related to sub-pixel variation.The root-mean-squared error (RMSE) for optical approaches is 0.08 to 0.16.Approaches using passive microwave data [19] are limited by the sensor resolution, especially in narrow channels where land contamination is likely.They are also only applicable to very high sea ice concentrations; even small areas of open water would dominate the signal.
Approaches using active microwave data have provided varying results.Barber and Yackel [20] report a relationship between ERS-1 backscatter and f p , but note a limitation related to wind roughening of melt ponds.Mäkynen et al. [21] report poor skill in connecting changes in f p with backscatter magnitude and texture using ENVISAT WSM images of sea ice during the ponding stages.Scharien et al. [22] report on the usefulness of co-polarization and cross-polarization ratios in estimating f p during ponding using RADARSAT-2 quad-pol data; their RMSE values range from 0.05 to 0.43.Using TerraSAR-X, Fors et al. [23] also find the co-polarization ratio useful during the ponding stage, at intermediate wind speeds, with an RMSE is 0.4; they also evaluate four polarimetric parameters and find reasonable correlations.To evaluate prediction of f p from winter images, Scharien et al. [24,25] use ENVISAT-ASAR and Sentinel-1 co-and cross-polarized backscatter and texture measures thereof.They report overall RMSE values of 0.08 and 0.09, but acknowledge model underestimation of f p for smooth FYI related to sensor noise floors.
Estimating melt pond fraction over smaller areas and for validation of the techniques above rely on in-situ observations [11,22], or low-level aerial photographs [26,27].A recent study using low-level aerial photographs showed that, due to their fine spatial and temporal resolutions, they are an ideal data source for estimating f p since image classification accuracies are as high as 97% [22].
Using winter imagery to predict f p implies an association of sea ice roughness and/or its snow cover variability with subsequent f p .Comprehensive studies on the degree of winter snow thickness spatial variability, as well as relationships with f p , are essential to better understand radiative transfer processes leading to snow and sea ice ablation rates during the spring melt transition.However, studying snow thickness variability on sea ice is complicated.The accumulation and redistribution of snow on FYI exhibits high spatio-temporal variability [9], and logistical challenges, resulting in a shortage of in-situ data [10].Snow thickness distributions have been obtained in-situ [9,[28][29][30], using laboratory-based procedures [31,32], and using remote sensing methods based on passive microwave data [33][34][35], frequency-modulated continuous-wave airborne radar returns [36,37] and Synthetic Aperture Radar (SAR) data [7,38,39], and using a combination of active microwave scatterometer and optical data [40].
The use of SAR backscatter must also contend with variability in backscatter with incidence angle (θ).Kim et al. [41] identified the θ effect on backscatter intensity, as it strongly influences surface and volume scattering processes.Microwaves are dominated by surface scattering at lower θ (<30 • ); whereas at higher θ (>30 • ), volume scattering dominates over surface scattering [42,43].
The advent of polarimetric SAR data provide additional parameters for characterizing sea ice surfaces.Polarimetric parameters such as the Co-polarized phase difference, Co-polarized correlation coefficient, Entropy, Anisotropy, and Alpha angle can be calculated from second-order derivatives of the scattering matrix (i.e., covariance and coherency) to obtain enhanced information about both surface and volume scattering mechanisms of snow and sea ice [39,[44][45][46][47].
In addition to linear and polarimetric parameters, SAR image texture also provides valuable structural information, leading to higher winter sea ice classification accuracy when compared to backscatter intensity alone [48][49][50][51].Second-order texture measures, derived from the gray-level co-occurrence matrix (GLCM) introduced by Haralick et al. [52], were evaluated for sea ice classification [49], showing that the combination of σ 0 (tone/grey-level) and texture measures (spatial dependence in tone) gives better results than using texture measures alone.GLCM texture is the most commonly used texture analysis technique; as it takes into account the spatial organization of grey tones within a moving window and offers a second-order statistical technique for extracting texture features [49,52].Image texture of snow covered sea ice is a function of its near surface characteristics e.g., snow properties, snow thickness variability, and ice roughness [53].Several studies successfully demonstrated the potential of GLCM texture measures, calculated from SAR images, to improve classification/segmentation of snow covered sea ice [49,51,53,54].
Given that the most consistent skill in estimating f p is reported by winter prediction of f p [24], we build on the winter prediction technique in this paper.We maintain the use of texture parameters as in Mäkynen et al. [21] and Scharien et al. [24] and expand on the analysis by evaluating the contributions of polarimetric parameters as per Fors et al. [23], and we assess texture measures of the polarimetric parameters.Furthermore, given the relationship of winter snow thickness and subsequent f p , we attempt an initial inversion of f p to estimate snow thickness.

Objectives
The primary objective of this study is to predict early summer f p from winter C-band SAR polarimetric backscatter.A secondary objective is to relate the predicted f p to the spatial distribution of snow thickness.Meeting these objectives will aid in understanding the relationship between snow thickness variability on FYI during winter and its evolution into a melt pond icescape in early summer.Since the SAR backscatter coefficient σ 0 is largely a function of θ, this study will also examine the θ dependency on backscatter and its relationship with f p .The following research questions towards these objectives are stated:

1.
What are the relationships between winter C-band RADARSAT-2 (RS-2) SAR backscatter (linear and polarimetric parameters and texture parameters) and f p ? 2.
How does the relationship between RS-2 SAR backscatter and f p change with θ? 3.
Can we predict f p based on linear, polarimetric and texture parameters, and can predicted f p be used to infer the late winter snow thickness variability?

Study Area
The study area is located in Resolute Passage, Nunavut, Canada and comprises of a large area of landfast FYI.The study area was imaged by the RADARSAT-2 SAR satellite during late-winter and by aerial photography during the melt pond season (Figure 1).Analyzing the melt pond evolution stages, we consider that stage 2 is a suitable time to examine aerial images of f p and its relationship with late winter snow thickness.In-situ observations were carried out from late-winter through the melt pond season at a field site ~4 km southwest of Sheringham Point, Cornwallis Island.Fieldwork was conducted in two coordinated sampling periods (Legs).Leg 1 field work was conducted during the late-winter/early-spring period (8 to 29 May 2012), while Leg 2 was during late-spring/summer (4 June to 2 July 2012).Melt-onset was observed on 1 June 2012.Throughout Leg 1, the snow thickness range in the study region was 1.5 cm to 40 cm.Snow depths from snow pits and 5 transect lines 100 to 150 m long were collected.One snow depth measurement was acquired from each pit (a total of 102), and 565 snow depth measurements were sampled from all transects at 1 m spatial intervals.For this study, we used 202 samples of snow depth, including 102 measurements from the snow pits, and 100 measurements randomly selected from the snow transects to minimize spatial autocorrelation (Figure 2a) [9].SAR acquisitions and snow samples were obtained at air temperatures between −20 • C and −3 • C (Figure 2b).An autonomous 2.75 m tall micro-meteorological station was installed on the sea ice from 18 May to 25 June 2012 to monitor air temperature, wind speed, and direction at 1-min intervals.
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 22 5 transect lines 100 to 150 m long were collected.One snow depth measurement was acquired from each pit (a total of 102), and 565 snow depth measurements were sampled from all transects at 1 m spatial intervals.For this study, we used 202 samples of snow depth, including 102 measurements from the snow pits, and 100 measurements randomly selected from the snow transects to minimize spatial autocorrelation (Figure 2a) [9].SAR acquisitions and snow samples were obtained at air temperatures between −20 °C and −3 °C (Figure 2b).An autonomous 2.75 m tall micro-meteorological station was installed on the sea ice from 18 May to 25 June 2012 to monitor air temperature, wind speed, and direction at 1-min intervals.

Aerial Surveys
In early summer, the field site was initially characterized by a very shallow and wet snow cover overlying a rough superimposed ice layer.Melt ponds started to form on 10 June 2012, were fully formed on 15 June reaching at a maximum pond fraction of 78%, and remained above 50% until 25 June [12].To quantify the evolution of melt ponds, their morphology, and f p at a regional scale, four aerial photography surveys were flown between 13 and 29 June 2012 via fixed-wing Twin Otter Remote Sens. 2018, 10, 1603 5 of 21 aircraft.A Canon G10 camera attached to an opening in the hindmost of the plane acquired digital optical images at nadir viewing angle, with a sampling rate resulting in 0% to ~10% image overlap.Data from the second survey flight on 22 June 2012, 01:43 UTC (Figure 1) are used in this study, while the melt ponds were fully formed and the geometric structures of the melt ponds were established.The survey was flown at an altitude of 610 m in a grid pattern that also coincided with the field study site.The image data has a swath width of 750 m and pixel size of 0.22 m [22].At the time of acquisition, the surface air temperature was approximately 2.6 • C, and the melt ponds were liquid and moderately wind roughened over the entire region; the wind speed was 3.2 m/s at the on-ice meteorological station and 5.3 m/s at the Resolute weather station.).This study uses seven RS-2 scenes acquired during Leg 1, that encompass the majority of aerial photographs of melt ponded sea ice acquired during Leg 2 (Table 1).2.3.Methods

Estimation of Melt Pond Fraction from Aerial Photographs
f p was calculated from aerial photographs using a decision-tree classification approach [22].Melt ponds are distinguished from adjacent ice/snow patches because the albedo of ponds is significantly lower than the snow and/or ice.The resulting classification (e.g., Figure 3a) has an error of 0.03 [22].The distribution of f p for the 22 June 2012 flightpath is presented in Figure 3b.

RADARSAT-2 Data Processing
RS-2 data were calibrated, and speckle filtered using a 5 × 5 boxcar filter to reduce the speckle effect while preserving image information [55,56].Four linear and five polarimetric parameters were derived (Table 2).Polarimetric parameters were derived from covariance (1) and coherency matrices (2), as applicable [57,58].Eight GLCM texture parameters (Table 3) were calculated for each linear and polarimetric parameter (Table 2) using a 5 × 5 moving window.Therefore, a total of 81 parameters (image bands) were generated.All parameter images were projected to a common geographic coordinate system using nearest neighbor intensity interpolation, at 5 m pixel spacing.

RADARSAT-2 Data Processing
RS-2 data were calibrated, and speckle filtered using a 5 × 5 boxcar filter to reduce the speckle effect while preserving image information [55,56].Four linear and five polarimetric parameters were derived (Table 2).Polarimetric parameters were derived from covariance (1) and coherency matrices (2), as applicable [57,58].Eight GLCM texture parameters (Table 3) were calculated for each linear and polarimetric parameter (Table 2) using a 5 × 5 moving window.Therefore, a total of 81 parameters (image bands) were generated.All parameter images were projected to a common geographic coordinate system using nearest neighbor intensity interpolation, at 5 m pixel spacing.

Image Sampling
Polygons corresponding to each 750 m by 750 m aerial photograph were used to extract spatially coincident SAR backscatter statistics from the seven RS-2 SAR scenes.The number of samples varied according to the portions of the RS-2 scenes that overlapped aerial photographs (Table 1).The total number of samples is 650, distributed across three different θ ranges.
Sample means of the nine linear and polarimetric parameters and eight GLCM parameters were calculated.Therefore, each of the 650 samples contains 81 associated mean values, corresponding to the nine linear and polarimetric parameters and their 72 derived texture parameters.Given that many adjacent aerial photographs overlap each other by up to approximately 10%, not all samples are entirely independent.However, the sample means and texture parameters represent the polygons as a whole, with no single-sample redundancy.Samples were divided into two equally sized training and test sets using a random sampling approach.A set of 325 training samples were used for model development, and 325 test samples were used to validate the derived models.

Terminology Theoretical Description References
Entropy (ENT) Detects the randomness of grey level distribution.Entropy will be higher for nearly random or noisy images.[59] Angular second moment (ASM) Detects disorder in texture.High energy values represent the constant of periodic form of gray level distribution [53,60] Contrast (CON) The difference between the highest and the lowest values of a connecting set of pixels.Contrast is greater for images that have rapidly fluctuating intensities. [61] GLCM variance (VAR) The degree of heterogeneity.Variance increases when the gray level values differ from their mean.[62] Correlation (COR) Measures linear-dependencies of gray tone in the image.[63] Homogeneity (HOM) Gives maximum value when all elements in the image are same.[49] Dissimilarity (DIS) Assists to distinguish most surface features.[64] GLCM Mean (MEAN) Measures the mean of Grey level co-occurrences values.Generally least correlated with other texture parameters. [65]

Regression Model Development
Multivariate linear regression was employed to generate combined backscatter, polarimetric parameter, and texture-based models to predict f p .Spearman's rank correlation analysis was used to assess correlation between parameters and with f p while accounting for any non-normal distributions.A coefficient threshold of 0.7 was used to identify multicollinearity between parameters.Correlated parameter pairs were assessed against f p , and the parameter with the higher correlation with f p was kept for regression model development.Model selection was then accomplished using forward stepwise regression analysis with an entry p-value of 0.05.These analyses were performed separately for the three θ ranges (NR, MR and FR).Therefore, a different model is derived for each of the three θ ranges.
To identify any misspecifications in the multivariate regression models, and to support the accuracy values obtained with multivariate regression, a random forest regression [66] was also performed, using the random Forest package in the statistical package-R.The number of parameters used at each split was set to the number of parameters divided by 3, and 500 trees were used.
The reduced parameter set, following removals due to correlation, was used for the random forest regression as well, since multicollinearity has been found to be detrimental [67].
Model accuracy was evaluated via RMSE analysis of predicted f p values versus observed values in the test set.This was done for both multivariate and random forest regressions.Multivariate regression models were used to estimate f p for RS-2 SAR scenes resampled to coarser 750 m by 750 m pixels to match the scale of the input sample data.To examine whether model estimates are sensitive to scale, images were also resampled to 500 m and 250 m pixel spacings.

Snow Thickness Images
Given the association of f p and winter snow thickness [7,8], we infer the spatial distribution of relative snow thickness based on predicted f p .To constrain the snow thickness, we assume that our in-situ-measured range of snow thickness (1.5 cm to 40 cm) is representative of the SAR scene extents (25 km by 25 km), which encompasses the study area.We linearly categorized the in-situ-measured snow thickness into eight classes, evenly-distributed across the snow thickness range, and related them to eight f p classes ranging from 0.25-0.85(i.e., lowest f p class is equivalent to the highest class of snow thickness).This provides a rudimentary hindcast of snow thickness.

Relationship between Winter SAR Backscatter and Pond Fraction
Spearman rank correlation coefficients between RS-2 parameters and f p are shown in Table 4.
HH is correlated with one or more parameters in each θ range and is thus not included in model development; it is shown in Table 4 for comparison.In the NR, σ 0 HV and MEAN R VV HH exceed 0.5, as does σ 0 HH .The value for σ 0 HV (−0.652) is within the range reported in Scharien et al. [24] (−0.647 to −0.792), but the value for σ 0 HH (−0.540) is slightly lower than their range (−0.590 to −0.824).Our relatively higher correlation value for σ 0 HV may relate to the lower noise floor of our data; suggesting that σ 0 HV was limited by the noise floor in Scharien et al. [24].Our study examines data at a much lower noise floor: −37 dB (NR) to −33 dB (FR) for RS-2 fine quad-polarization mode (MDA, 2014), versus −20 dB to −22 dB for Envisat-ASAR and Sentinel-1 [68,69].For σ 0 HV , our mean backscatter exceeds the RS-2 fine quad-polarization noise floor by 6.4 dB at NR, by 4.0 dB at MR, and by 3.1 dB at FR.In the MR, three parameters exceed 0.5: MEAN HV , MEAN HH and ASM HV , plus σ 0 HH .Our correlation for σ 0 HH (−0.552) is lower than the range reported by Scharien et al. [24] (−0.695 to −0.810).However, our lower noise floor is likely the cause of our texture parameters derived from σ 0 HV exhibiting relatively strong correlation with f p at MR when using RS-2, whereas the noise floor was found to be a limiting factor in Scharien et al. [24].In the FR, three parameters exceed 0.5: σ 0 VV , VAR VV and MEAN A , plus σ 0 HH .We find σ 0 VV to consistently have higher correlation with f p than σ 0 HH , suggesting that employing σ 0 VV instead of σ 0 HH may improve the results of Scharien et al. [24].Overall, our correlation values are not as high as those found by Scharien et al. [24], which may be the result of their object-based approach versus our polygon-based approach.In comparison to Fors et al. [23], we find higher correlations for H, α, ρ and φ at FR, but lower correlations for H, α, ρ and φ at NR.Given that their X-band correlations are for sea ice with existing melt ponds, versus our C-band correlations for snow covered sea ice during winter, discrepancies are not unexpected.
The correlation analysis reveals multicollinearity between many parameters (r ≥ 0.7).Of the original 81 parameters, 64 are correlated at NR, 57 at MR and 62 at FR; the remaining less-correlated parameters are used for model development; these are shown in Table 4.
Given that late winter σ 0 VV and σ 0 HV exhibit negative correlation with melt season f p (Table 4; e.g., Figures 4 and 5), and that f p is inversely correlated with late winter snow thickness [7], we infer that σ 0 VV , and σ 0 HV increase with snow thickness in a manner similar to previous studies [2,39].R VV HH exhibits substantial negative correlation with f p at MR (e.g., Figure 6).R VV HH pixel values are primarily negative, indicating that σ 0 HH > σ 0 VV , suggesting that Fresnel reflection effects associated with a smooth dielectric interface are responsible [70].R VV HH tends toward 0 dB as f p decreases, pointing to reduced Fresnel reflection, which, in turn, may be associated with rougher FYI and thicker snow.with  are sensitive to ice and snow surface structure and related melt season  .Texture parameters based on  and  are present in all three  ranges, suggesting that these parameters may be sensitive to snow thickness differences and are therefore correlated with  .Texture parameters most correlated with f p are primarily based on σ 0 HV , σ 0 VV σ 0 HH , R VV HH , φ and A (Table 4, Figures 7-9).Texture parameters that are correlated with f p and derived from σ 0 HV are present in all three θ ranges.Their correlation with f p suggests that σ o HV might be sensitive to different snow thicknesses as well.σ o HV is influenced by surface roughness, volume scattering and multi-bounce scattering [71]; as snow thickness increases, the potential for increased volume scattering increases [72].Thicker snow covers are associated with rougher sea ice [9], which may also influence σ o HV because cross polarization is the result of multiple volume scattering enhanced by the presence of uneven surfaces [73].Therefore, it is not surprising that texture parameters associated with σ 0 HV are sensitive to ice and snow surface structure and related melt season f p .Texture parameters based on φ and A are present in all three θ ranges, suggesting that these parameters may be sensitive to snow thickness differences and are therefore correlated with f p .

Prediction of Pond Fraction
The regression analysis produced the following models for predicting  at NR, MR and FR, respectively:

Prediction of Pond Fraction
The regression analysis produced the following models for predicting  at NR, MR and FR, respectively:

Prediction of Pond Fraction
The regression analysis produced the following models for predicting f p at NR, MR and FR, respectively: The models, as represented by Equations ( 3)-( 5), exhibit R 2 values of 0.57, 0.61, and 0.62 respectively.The FR model exhibits the highest R 2 and contains a linear parameter and four texture parameters.The MR model exhibits the second highest R 2 and contains one linear parameter and four texture parameters.All models make use of polarimetric data for the contributing texture parameters.The robustness of the models is consistent with previous findings in that the inclusion of texture with linear and/or polarimetric parameters improves classification results [48][49][50][51].
Differing scattering mechanisms at NR, MR and FR are very likely responsible for the three different models.Surface roughness changes depending on the local θ [74]; surface scattering dominates at NR, and volume scattering dominates at FR [41,43].In addition, the polarimetric parameters that were found to have low correlation with other parameters are likely sensitive to specific scattering mechanisms; this may explain their inclusion in the models even if their correlation with f p is low (Table 4).In comparison to the models developed by Scharien et al. [24], the model parameters here suggest that the addition of polarimetric capabilities can be a benefit, and a lower noise floor is a benefit when using HV.However, the use of additional model parameters may lead to over-fitting.
Random forest regression analysis provides predictive ability similar to the linear regression (see Section 3.3).An assessment of random forest parameter importance supports the parameter selection in the linear regression models (Table 5); most of the linear regression parameters are found in the top five important positions.Differences between models are expected due to their different selection techniques.Give the random forest results, no gross misspecification is anticipated in the linear models.

Error Analysis
Predicted f p , using Equations ( 3)-( 5), is compared with measured f p for the 325 test samples, for each θ range (Figure 10).The RMSE was found to be 0.106 for NR, 0.108 for MR, and 0.094 for FR.
The RMSE values obtained with random forest regression are 0.109, 0.091 and 0.092, for NR, MR and FR models, respectively.The NR linear model exhibits greater skill than the random forest model; however, the MR linear model is substantially worse that for random forest.Nevertheless, these values indicate that the multivariate linear regression models have predictive capabilities similar to that of random forest models and can thus be used without significant loss of skill.

SAR Scene Pond Fraction Estimates
The NR, MR and FR  prediction models are applied to late winter SAR scenes with pixel spacing of 750 m, 500 m and 250 m to assess the dependency of pixel size (Figures 11-13).The areas represented by the models at lower and higher resolution indicate similar pattern and variability, indicating that the models are appropriate at higher resolution.The higher resolution SAR scene displays higher variability than the 750 m or 500 m scenes, which is logical as 250 m pixel spacing would reveal more detail.There are possible mixed-pixel effects where sea ice abuts land because the land mask used may not be perfectly accurate in all cases.However, these effects only impact the sea ice pixels immediately adjacent to the coastlines.Buffering the coastline would mitigate this issue, if needed.
Within the area of overlapping scenes (red outline in Figure 11a), some differences in the estimates can be observed.In general, the NR estimates (Figure 11) exhibit lower pond fractions than do the MR (Figure 12) and FR (Figure 13) estimates; especially at the north end of Griffith Island and associated with the linear feature extending from Sheringham Point.The latter is a region of rougher, ridged, sea ice which is expected to have thicker snow cover and smaller  .The efficacy of the NR model to show such smaller  is supported by the error analysis which shows that NR exhibits somewhat less bias at small  (Figure 10a).Conversely, the FR estimates show larger areas of higher  than MR and NR.This is also supported by the error analysis, which shows a smaller bias at high  for FR (Figure 10c).The MR estimates appear to lack the range of  , compared to the NR and FR estimates; this is also displayed in the error analysis (Figure 10b).Given these observations, it may be beneficial to make use of combined NR and FR estimates.The estimates of  with regard to incidence angle are very likely a function of the different scattering mechanisms that dominate at the respective incidence angle ranges and the sensitivity of the different parameter combinations in the models to these scattering mechanisms.It is apparent that the models overestimate lower f p and underestimate higher f p Figure 10.The underestimation of higher f p is most pronounced at MR (Figure 10b).However, the linear dependencies are consistent across the three θ ranges.

SAR Scene Pond Fraction Estimates
The NR, MR and FR f p prediction models are applied to late winter SAR scenes with pixel spacing of 750 m, 500 m and 250 m to assess the dependency of pixel size (Figures 11-13).The areas represented by the models at lower and higher resolution indicate similar pattern and variability, indicating that the models are appropriate at higher resolution.The higher resolution SAR scene displays higher variability than the 750 m or 500 m scenes, which is logical as 250 m pixel spacing would reveal more detail.There are possible mixed-pixel effects where sea ice abuts land because the land mask used may not be perfectly accurate in all cases.However, these effects only impact the sea ice pixels immediately adjacent to the coastlines.Buffering the coastline would mitigate this issue, if needed.
Within the area of overlapping scenes (red outline in Figure 11a), some differences in the estimates can be observed.In general, the NR estimates (Figure 11) exhibit lower pond fractions than do the MR (Figure 12) and FR (Figure 13) estimates; especially at the north end of Griffith Island and associated with the linear feature extending from Sheringham Point.The latter is a region of rougher, ridged, sea ice which is expected to have thicker snow cover and smaller f p .The efficacy of the NR model to show such smaller f p is supported by the error analysis which shows that NR exhibits somewhat less bias at small f p (Figure 10a).Conversely, the FR estimates show larger areas of higher f p than MR and NR.This is also supported by the error analysis, which shows a smaller bias at high f p for FR (Figure 10c).The MR estimates appear to lack the range of f p , compared to the NR and FR estimates; this is also displayed in the error analysis (Figure 10b).Given these observations, it may be beneficial to make use of combined NR and FR estimates.The estimates of f p with regard to incidence angle are very likely a function of the different scattering mechanisms that dominate at the respective incidence angle ranges and the sensitivity of the different parameter combinations in the models to these scattering mechanisms.

SAR-Based Estimation of Snow Thickness
The snow thickness of our in-situ sample site ranges from 1.5 cm to 40 cm, with a mean thickness of 10.9 cm ± 7.7 cm (Figure 2).This range of snow thickness is representative of FYI in the Canadian Arctic Archipelago; however, the mean snow thickness is lower than values previously found for relatively smooth FYI [75].Considering the range of snow thickness in our study area, first-order hindcasts of snow thickness distribution were generated (Figures 14-16

SAR-Based Estimation of Snow Thickness
The snow thickness of our in-situ sample site ranges from 1.5 cm to 40 cm, with a mean thickness of 10.9 cm ± 7.7 cm (Figure 2).This range of snow thickness is representative of FYI in the Canadian Arctic Archipelago; however, the mean snow thickness is lower than values previously found for relatively smooth FYI [75].Considering the range of snow thickness in our study area, first-order hindcasts of snow thickness distribution were generated (Figures 14-16

SAR-Based Estimation of Snow Thickness
The snow thickness of our in-situ sample site ranges from 1.5 cm to 40 cm, with a mean thickness of 10.9 cm ± 7.7 cm (Figure 2).This range of snow thickness is representative of FYI in the Canadian Arctic Archipelago; however, the mean snow thickness is lower than values previously found for relatively smooth FYI [75].Considering the range of snow thickness in our study area, first-order hindcasts of snow thickness distribution were generated (Figures 14-16) based on the association of snow thickness and f p .The prediction maps provide insight into the snow thickness distribution within a 25 km by 25 km SAR scene.Of note is the higher accumulation of snow near the coastlines.This is because during wind drift, the uneven land surface topography causes more snow to accumulate on the leeward side [76] than in the homogeneous and smooth areas of FYI.The models illustrate a similar pattern of snow thickness distribution and show snow thickness variability at both fine and coarser resolution.Snow accumulation also depends on ice deformation [77], with snow catchment higher in deformed ice areas compared to smooth ice surfaces.The model hindcast maps indicating appreciably higher snow thickness along linear features near to the Cornwallis Island coastline appear to be reasonably well associated with the spatial location of narrow sea ice ridges as observed in the winter SAR imagery.Most of the study area exhibits a thin snow cover consistent with a relatively smooth surface.At both 750 m and 250 m pixel spacing, the models show similar pattern and variability of snow thickness.within a 25 km by 25 km SAR scene.Of note is the higher accumulation of snow near the coastlines.This is because during wind drift, the uneven land surface topography causes more snow to accumulate on the leeward side [76] than in the homogeneous and smooth areas of FYI.The models illustrate a similar pattern of snow thickness distribution and show snow thickness variability at both fine and coarser resolution.Snow accumulation also depends on ice deformation [77], with snow catchment higher in deformed ice areas compared to smooth ice surfaces.The model hindcast maps indicating appreciably higher snow thickness along linear features near to the Cornwallis Island coastline appear to be reasonably well associated with the spatial location of narrow sea ice ridges as observed in the winter SAR imagery.Most of the study area exhibits a thin snow cover consistent with a relatively smooth surface.At both 750 m and 250 m pixel spacing, the models show similar pattern and variability of snow thickness.
Between the NR, MR and FR snow thickness estimates, NR (Figure 14) exhibits the greatest snow depths.This relates directly to the smaller  estimates with the NR model as discussed in the previous section.The reduced bias at small  with NR (Figure 10a) likely provides the best estimates of thick snow.Thinner snow estimates appear to have the most realistic spatial distribution at FR (Figure 16).Again, this relates directly to the smaller bias at large  with FR (Figure 10c), where thinner snow is expected.

Conclusions
This paper presents an approach utilizing late winter C-band RADARSAT-2 SAR linear, polarimetric and texture parameters of snow covered first-year sea ice to predict early summer melt pond fraction.Predictions were then used to examine relationships between early melt pond fraction and late winter snow thickness.We estimated melt pond fractions from the aerial photographs collected over first-year sea ice.The correlation of melt pond fraction and late winter backscatter parameters provided us with to the ability to hindcast snow thickness.Nine linear and polarimetric parameters and 72 texture parameters were evaluated for their relationship with pond fraction.
Multivariate models comprised of linear, polarimetric and/or texture parameters are derived at near-, mid-and far-range incidence angles.The best pond fraction prediction capability is exhibited by the model at far-range incidence angles (RMSE = 0.09).Between the NR, MR and FR snow thickness estimates, NR (Figure 14) exhibits the greatest snow depths.This relates directly to the smaller f p estimates with the NR model as discussed in the previous section.The reduced bias at small f p with NR (Figure 10a) likely provides the best estimates of thick snow.Thinner snow estimates appear to have the most realistic spatial distribution at FR (Figure 16).Again, this relates directly to the smaller bias at large f p with FR (Figure 10c), where thinner snow is expected.

Conclusions
This paper presents an approach utilizing late winter C-band RADARSAT-2 SAR linear, polarimetric and texture parameters of snow covered first-year sea ice to predict early summer melt pond fraction.Predictions were then used to examine relationships between early melt pond fraction and late winter snow thickness.We estimated melt pond fractions from the aerial photographs collected over first-year sea ice.The correlation of melt pond fraction and late winter backscatter parameters provided us with to the ability to hindcast snow thickness.Nine linear and polarimetric parameters and 72 texture parameters were evaluated for their relationship with pond fraction.
Multivariate models comprised of linear, polarimetric and/or texture parameters are derived at near-, mid-and far-range incidence angles.The best pond fraction prediction capability is exhibited by the model at far-range incidence angles (RMSE = 0.09).
By relating late winter snow thickness and f p , these predictions help us to understand the winter snow thickness variability.The models are able to distinguish higher snow thickness along sea ice ridges, coastlines and relatively thinner snow cover on smooth surfaces of first-year sea ice, which is consistent with previous findings.
Moreover, the models show that the combination of SAR linear and polarimetric backscatter and texture parameters enhance the strength of the models compared to utilizing them separately for the prediction of pond fraction.The estimation of pond fraction over an entire SAR scene based on the models show logical distributions of melt ponds and snow thickness.
The results of this study add to the suite of seasonal sea ice forecasting tools, and thus can aid ship navigability since melt ponds are associated with advanced sea ice melt and significant weakening of sea ice mechanical strength.At the same time, it provides insight into the late winter snow thickness distribution on first year sea ice.The method is tested on landfast sea ice, is sensitive to the time periods of the collected aerial melt pond distributions and SAR scene acquisitions, and is likely not applicable to drift ice.Future research should test this model on a regional scale, and similar models should also be evaluated for their application over multi-year sea ice and more deformed types of first-year sea ice.
Author Contributions: S.R., T.G. and J.Y. designed the experiment.S.R. formulated the research methodology and wrote the manuscript.R.S. provided necessary data and contributed to manuscript development and proof reading.T.G. and J.Y. contributed with additional inputs during manuscript development and proof reading.

Figure 1 .Figure 1 .
Figure 1.Study area in Resolute Passage, Nunavut, illustrating the in-situ field site, the flightpath of the aerial survey, and the region covered by RADARSAT-2 scenes.The red box in the inset map of the Canadian Arctic Archipelago delineates the location of the study area.

Figure 1 .
Figure 1.Study area in Resolute Passage, Nunavut, illustrating the in-situ field site, the flightpath of the aerial survey, and the region covered by RADARSAT-2 scenes.The red box in the inset map of the Canadian Arctic Archipelago delineates the location of the study area.

Figure 2 .
Figure 2. (a) Distribution of 202 in-situ snow thickness measurements collected on FYI acquired from the field site during Leg 1.(b) Hourly air temperature at the Resolute CARS weather station, Nunavut (blue), and once a minute on-ice air temperature (green).2.2.2.RADARSAT-2 Data RS-2 SAR scenes were acquired in fine quad-pol mode over the aerial flightpath at near-range (18 • -29 • ) (NR), mid-range (30 • -39 • ) (MR) and far-range (40 • -49 • ) (FR) θ.Each scene consists of a polarimetric (HH + VV + HV + VH and inter-channel phase information) data set, with a nominal 5.2 m × 7.7 m resolution in range and azimuth, respectively, and covers a 25 km × 25 km area (MDA, 2014).This study uses seven RS-2 scenes acquired during Leg 1, that encompass the majority of aerial photographs of melt ponded sea ice acquired during Leg 2 (Table1).

Figure 3 .
Figure 3. (a) Sample of a classified aerial photograph with melt ponds (grey areas) and snow patches (white areas); (b) histogram of pond fraction distribution for the flightpath on 22 June 2012.

Figure 3 .
Figure 3. (a) Sample of a classified aerial photograph with melt ponds (grey areas) and snow patches (white areas); (b) histogram of pond fraction distribution for the flightpath on 22 June 2012.

Figure 4 .
Figure 4. Late winter  dB as a function of  at NR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 4 . 22 Figure 5 .
Figure 4. Late winter σ 0HV (dB) as a function of f p at NR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 6 .
Figure 6.Late winter  dB as a function of  at MR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 7 . 63 Figure 5 . 22 Figure 5 .
Figure 7. Late winter  , as a function of  at MR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 6 .
Figure 6.Late winter  dB as a function of  at MR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 7 . 63 Figure 6 .
Figure 7. Late winter  , as a function of  at MR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 6 .
Figure 6.Late winter  dB as a function of  at MR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 7 . 63 Figure 7 . 22 Figure 8 .
Figure 7. Late winter  , as a function of  at MR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 9 .
Figure 9. Late winter  , as a function of  at NR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 8 . 22 Figure 8 .
Figure 8.Late winter MEAN φ , as a function of f p at FR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 9 .
Figure 9. Late winter  , as a function of  at NR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 9 .
Figure 9. Late winter COR A , as a function of f p at NR.The solid line represents the linear dependency, and r is the Spearman correlation.

Figure 10 .
Figure 10.Observed (estimated via aerial photography) and regression model predicted  at (a) NR, (b) MR and (c) FR.The dashed line is the linear dependency of the comparison.The solid line represents the 1:1 relationship.

Figure 10 .
Figure 10.Observed (estimated via aerial photography) and regression model predicted f p at (a) NR, (b) MR and (c) FR.The dashed line is the linear dependency of the comparison.The solid line represents the 1:1 relationship.

Figure 15 .
Figure 15.Snow thickness class spatial distribution from the MR model at 750 m (a), and 250 m (b) pixel spacing, respectively.Grey area represents land.Based on RS-2 FQ15 2012-05-21 data.

Figure 16 .
Figure 16.Snow thickness class spatial distribution from the FR model at 750 m (a) and 250 m (b) pixel spacing.Grey area represents land.Based on RS-2 FQ23 2012-05-15 data.

Figure 16 .
Figure 16.Snow thickness class spatial distribution from the FR model at 750 m (a) and 250 m (b) pixel spacing.Grey area represents land.Based on RS-2 FQ23 2012-05-15 data.

Funding:
This research was funded by Canadian NSERC Discovery grants to John Yackel and Randy Scharien as well as MEOPAR (Marine Environmental Observation Prediction and Response Network) funding to Randy Scharien.The APC was funded by Canadian NSERC Discovery grants to John Yackel.

Table 1 .
25km by 25 km RS-2 scenes used in this study and the number of 750 m by 750 m aerial photographs that fall within each SAR scene.

Table 3 .
GLCM texture parameters used in the study.

Table 4 .
Parameters used for model development.Parameters have low correlation with each other (Spearman's r < 0.7).Parameters are shown in descending order of absolute correlation (Spearman's) with f p .NR = near-range, MR = mid-range, and FR = far-range incidence angles.σ 0 HH (value in italics) is included for comparison.

Table 4 .
Parameters used for model development.Parameters have low correlation with each other (Spearman's r < 0.7).Parameters are shown in descending order of absolute correlation (Spearman's) with  .NR = near-range, MR = mid-range, and FR = far-range incidence angles. (value in italics) is included for comparison.

Table 5 .
Random forest regression parameters in order of parameter importance (increase in mean-squared-error with removal (IncMSE)), for NR, MR and FR models.The top 13 important positions are shown.Shaded parameters are found in the respective linear regression models.