Tenuous Correlation between Snow Depth or Sea Ice Thickness and C- or X-Band Backscattering in Nunavik Fjords of the Hudson Strait

: Radar penetration in brine-wetted snow-covered sea ice is almost nil, yet reports exist of a correlation between snow depth or ice thickness and SAR parameters. This article presents a description of snow depth and ﬁrst-year sea ice thickness distributions in three fjords of the Hudson Strait and of their tenuous correlation with SAR backscattering in the C- and X-band. Snow depth and ice thickness were directly measured in three fjords of the Hudson Strait from 2015 to 2018 in April or May. Bayesian linear regression analysis was used to investigate their relationship with RADARSAT-2 (C-band) or TerraSAR-X (X-band). Polarimetric ratios and the Cloude–Pottier decomposition parameters were explored along with the HH, HV and VV bands. Linear correlations were generally no higher than 0.3 except for a special case in May 2018. The co-polarization ratio did not perform better than the backscattering coefﬁcients.


Context
Seasonal snow and ice covers in Nunavik are affected by the impacts of climate change: Kangiqsujuammiut (people of Kangiqsujuaq, Nunavik, in Canada) have reported later sea ice freeze-up in the fall [1], as well as less snow on the ground, earlier sea ice breakup in spring, changes in travel routes and more variable winds [2]. Inuit lives are embedded in the climate change context, and its impacts on sea ice have practical and immediate consequences on personal safety and access to travel and marine wildlife [3]. With shipping traffic in the Canadian Arctic having markedly increased over the last decade [4] and 2040-2064 climate projections for the region showing shorter snow cover periods [5], sea ice conditions and their impact on land-use and marine transport will continue to evolve. Yet, historically, the scientific community has undertaken monitoring efforts at scales too coarse to account for regional or local variations in ice conditions [3].
In this study, we investigate snow depth and ice thickness distributions in Salluit, Deception Bay, and Kangiqsujuaq, all Nunavik fjords of the Hudson Strait, as well as their correlation with C-and X-band SAR. This work is relevant due to land use and shippingrelated operations by communities and industries in the area and for the sea ice remote sensing community. It is part of the Ice Monitoring project, a collaboration between the

Relationship with C-and X-Band Backscattering
SAR sensors are widely used for sea ice monitoring, both by governments and in scientific communities, due to their ability to operate in the presence of clouds and at night-time [11]. Current applications of this technology include characterizing melt pond fraction [12,13] and community-relevant sea ice features [10,14], as well as estimating ice stability [15]. Recent work by Yackel et al. [16] showed that, prior to melting onset, the daily variance in C-and Ku-band backscattering is correlated with relative snow thickness. While the majority of these sea ice applications rely on the C-band [11], other SAR frequencies like the Ku-, X-, and L-bands have also been investigated [17][18][19][20]. The X-band can adequately discriminate between newly formed ice and its surroundings [17] and is more sensitive to melt onset and surface roughness than the C-band [21], as well as changes in top-layer snow salinity [19].
This article and the references cited hereafter focus on undeformed first-year sea ice in cold conditions with air temperatures below the freezing degree point. We explore the VV, HV and VV polarizations as well as polarimetric ratios and the Cloude-Pottier decomposition parameters [22]. Nandan et al. [23] have recently demonstrated the extent to which SAR penetration in the snow covering first-year sea ice is limited by brine. They reported a maximum penetration of four cm into the brine-wetted snow for the C-and X-bands, near Resolute Bay in Nunavut. How can these results be reconciled with several observations of a correlation between undeformed sea ice thickness or snow depth and SAR parameters? Indeed, linear correlation coefficients of 0.4 to 0.6 were reported between ice thicknesses ranging from 15 to 150 cm and the HH or VV backscattering coefficients in the X-, C-or L-band over undeformed first-year drift ice in the Sea of Okhotsk [24,25]. In these cases, the co-polarization ratio (VV/HH) performed better than the individual polarizations, reaching an r-squared of 0.8 for the C-band [25]. A strong exponential relationship was observed between ice thickness ranging from 0 to 200 cm and the C-band CP ratio for drift and landfast undeformed first-year ice in the Labrador Sea, based on simulated compact polarimetry data [26]. This ice thickness estimation parameter is based on circular right-hand transmission and linear HH and VV detection, i.e., the compact polarimetry mode developped for the RADARSAT Constellation Mission [26].
As for snow covers on smooth landfast first-year ice, Gill et al. [27] noted that a positive correlation between C-band HH backscattering and snow depth is best observed in cold (−7.9 • C) than warm conditions (−0.4 • C). The associated r-squared values ranged from 0.3 to 0.8 for incidence angles of 27 • to 36 • , with measurements performed in Franklin Bay, Northwest Territories. In contrast, Nandan et al. [19] rather observed that the HH backscattering was greater from a thin snow cover of 4 cm than from snow 8 or 14 cm deep, which they attributed to steeper salinity gradients in thin snow causing enhanced surface scattering. In this study of sea ice near Resolute Bay, Nunavut, they also noted that the HH polarization was more sensitive to snow depth variations than its VV counterpart.

Objectives
This article combines field measurements and remote sensing data in order to characterize seasonal snow-covered sea ice from 2015 to 2018 in three Nunavik fjords of the Hudson Strait. The objectives of this article are to (1) characterize and explain snow depth and ice thickness distributions in Salluit, Deception Bay, and Kangiqsujuaq over three winters (2015-2018) and (2) investigate the empirical relationship between C-and X-band polarimetric SAR parameters and snow depth and ice thickness.

Study Areas
The study sites are three neighbouring fjords located along the coast of Nunavik in the Hudson Strait ( Figure 1). Two of them are home to Inuit communities: Salluit and Kangiqsujuaq. The third, Deception Bay, is located between the two communities, 50 km west of Salluit. Two mining companies have the marine infrastructure in Deception Bay, where their icebreakers transit year-round except from the mid-March to 1 June black-out window [28].
Freeze-up occurs in November or December and breakup in June or July [1,[29][30][31]; freeze-up generally proceeds through consolidation of young ice types. During the ice season, the study areas are covered in undeformed and smooth landfast first-year sea ice, except along the broken ice track left by ice-breaking transport in Deception Bay and near the shores where the tides lead to deformation. These features were excluded from the study. There is typically some small-scale surface roughness, on the order of the radar wavelength. Some areas may also sometimes feature either radar-smooth ice, or ice pieces sticking out of the ice by up to 15 cm. Monthly total precipitation data from Salluit airport is presented in the supplementary materials ( Figure S1). No weather station was operating in Deception Bay during the study and no precipitation data was available from the Kangiqsujuaq airport weather station. The bathymetry maps ( Figure S2 in the supplementary) show maximum water depths greater than 100 m in Salluit, 80 m in Deception Bay, and 200 m in Kangiqsujuaq. The difference between high and low tide ranges from 1.5 to 5 m in Salluit, from 2 to 5.5 m in Deception Bay, and from 4 to 8.5 m in Kangiqsujuaq [32]. In Deception Bay, GENIVAR [28] measured water salinity between 29 and 33 ppt. Top-layer ice salinity measured in January 2018 ranged from 5 to 10 psu [29]. Brine-wetted snow salinity near the ice surface ranged from 23 to 33 psu in January 2018 [29] and from 7 to 19 psu in May 2018.

Snow Depth and Sea Ice Thickness Measurements
Measurements were performed in January-February and April-May of 2016, 2017 and 2018, for each site, except for April 2016 where bad weather prevailed in Salluit (Table 1). Snow depth was measured using a meter-rule and ice thickness using a 2.5 cm diameter Kovacs ice auger and measuring tape. For a given site, sampling was done at 20 to 30 target locations arranged in a grid-like pattern with 1 to 2 km spacing ( Figure 2) and within 48 h or less except for Kangiqsujuaq in May 2018. This article focuses on the end-of-winter April-May data.    Tables 2 and 3 show mean snow depth and ice thickness for January-February and April-May measurements, respectively. Kangiqsujuaq always presented the deepest snow and thinnest ice (except in January 2017). Snow was deeper in January 2016 than in 2018 (Table 2), despite earlier measurements by one week (Table 1). Ice thickness at the end of the winter (Table 3) was greater in 2016 than in the other two years for both Deception Bay and Kangiqsujuaq; no data are available for Salluit in April-May 2016. Deception Bay presented its thinnest ice cover in 2017, whereas both 2017 and 2018 were similar for Salluit and Kangiqsujuaq.
Note that even in 2016 with an average thickness of 1.43 m, the ice in Deception Bay fell short of the historical range. Indeed, the most recent measurements, probably dating back to 1991, gave a thickness of 1.7 to 2 m [33]. In contrast, thicknesses measured from 2015 to 2018 ranged from 1.10 to 1.65 for extreme values and 1.20 to 1.43 on average. The ice in Deception Bay was therefore generally 50 cm thinner during our study than 25 years prior. Although historical data for the other two sites was not available, reports from 2007 by local experts stated that the ice was forming later and growing thicker more slowly than before in Salluit and Deception Bay [33].  The Canadian Ice Service secured RADARSAT-2 Wide-Fine Quad-Pol single look complex images for each site from 2015 to 2018 (Table 4). The C-band satellite operates at 5.405 GHz (5.55 cm wavelength), with a repeat period of 24 days. The acquisitions were performed in the descending orbit at roughly 6:30 local time (UTC-5 h). Six images were acquired per winter except in 2016, where a Deception Bay acquisition was cancelled. Scene size was 50 by 25 km before subsetting, with a spatial resolution of 5.2 and 7.6 m (range and azimuth, respectively) [34]. Air temperature at the time of acquisition was always below −10 • C except for Kangiqsujuaq in May 2018 where it was −3 • C (see Tables S1 and S2 in the supplementary). The noise-equivalent σ 0 value in the Wide-Fine Quad-Pol mode is −33± 6 dB [34].
The RADARSAT-2 images were processed with SNAP's Sentinel-1 Toolbox (version 6.0.0), a European Space Agency (ESA) open-access software. The toolbox was used through java-snap, a Java application made available on Gitlab [35]. After subsetting to the study areas, the data was converted from digital number to backscattering coefficient σ 0 [36]. Speckle filtering was performed over a 7 × 7 window using the Refined Lee or Polarimetric Refined Lee filter depending on the desired output (backscattering coefficient or covariance matrix). Lee filters preserve image details and contours [37]. This processing follows the ESA Polarimetric Tutorial [38], in which no multilooking is performed. We computed the Cloude-Pottier decomposition to extract the entropy (H), anisotropy (A) and alpha angle (alpha) [22,39]. Data position in the H-alpha plane provides information on the nature of dominant scattering mechanisms, for instance, surface or volume scattering [39]. The geometric correction was then performed using the Range-Doppler Terrain Correction algorithm using the freely accessible Canadian Digital Elevation Model [40] and nearest neighbour resampling for both the image and the DEM. Pixel spacing after this step was 8 m. The images did not present a large enough homogeneous area to estimate the equivalent number of looks after processing using the procedure outlined by Anfinsen et al. [41]. After geometric correction, linear σ 0 was converted to decibels. The co-polarization VV/HH and the cross-polarization VH/VV and HV/HH ratios were computed from linear σ 0 .

TerraSAR-X
The DLR secured TerraSAR-X or TanDEM-X acquisitions of StripMap dual-pol single look complex images over Deception Bay from 2015 to 2018. Acquired in ascending orbit 13 at 17:30 local time (UTC-5 h) with polarizations HH and VV, the images have an incidence angle of 38 • . A total of 75 images were acquired between 2015-12-23 and 2018-07-26. The X-band satellites operate at 9.65 GHz (3.11 cm wavelength) with a repeat period of 11 days. Images were acquired in a dual-pol orbit with polarizations HH and VV. The scene size before subsetting to the study area was 15 by 50 km, with a spatial resolution of 0.9 and 2.5 m for range and azimuth, respectively [42]. Air temperature at the time of acquisition was always below −10 • C (see Table S2 in the supplementary). The noise-equivalent σ 0 value in the StripMap mode for orbit 13 is −24.5 dB [42].
The DLR data was processed by using their in-house Multi-SAR System [43]. This processing includes converting from digital number to backscattering coefficient, multilooking to produce square pixels and increase radiometric quality, geometric correction using bilinear interpolation for the DEM and cubic convolution resampling for the image, and image enhancement to reduce speckle [44]. The output images have a pixel spacing of 2.5 m pixels with a radiometric resolution of 1.6 looks. Linear σ 0 was converted to decibels.

Computing SAR Parameter Seasonal Medians
Seasonal medians were computed from the satellite images acquired in each season. Areas of interest (AOIs) roughly 120 by 100 m and each containing between 600 and 650 pixels were distributed over the homogeneous study areas in a grid-like pattern with 0.7 to 1 km separation, avoiding special features like the shore or a ship's track. The Salluit, Deception Bay, and Kangiqsujuaq study areas counted 35, 43, and 78 AOIs, respectively. Median backscattering was computed over each AOI and then over all AOIs for a given image, yielding a single median value per image. Only images between January and May were used to avoid the freeze-up and spring transitions. The seasonal medians were computed from the resulting time-series to yield a single value per year, allowing interannual comparisons between sites. This step was performed using Python [45]. For the comparison with the X-band, this step was also performed on the processed TerraSAR-X images as described in Dufour-Beauséjour et al [31]; 32 of the 43 AOIs were covered by the TerraSAR-X images, with a higher spatial resolution.

Extracting SAR Parameter Values Coincident with Thickness Measurements
For every field campaign except those in January 2016, thickness measurements were paired with the image acquired in closest temporal proximity. The relative timing of fieldwork and image acquisitions is presented in the supplementary materials (Tables S1  and S2). We used Python to identify pixels coinciding with sampling locations and to extract the backscattering coefficient value over that location [46]. To mitigate speckle noise, a 3 × 3 pixel window average was used for the RADARSAT-2 images, corresponding to an area size of 576 m 2 . For the comparison between the C-band and the X-band, mean TerraSAR-X values were extracted over a 9 × 9 pixel window, which corresponds to 506 m 2 . The co-polarized SAR backscattering values were well above the noise-equivalent σ 0 for each sensor. Linear correlation coefficients were computed between thickness variables and all available SAR parameters. Because the Bayesian linear regression is a relatively new method in our field, we chose to apply it to a single parameter. The HH backscattering coefficient was chosen because it is commonly used in sea ice applications, was available for both radar frequencies, and was not systematically outperformed by another parameter in terms of the linear correlation coefficient (see Table S3 in the supplementary).

Bayesian Framework
We used Bayesian linear regression to investigate a potential linear relationship between snow depth or ice thickness and C-or X-band SAR. The advantage of the Bayesian framework, compared to the frenquentist approach (e.g., Pearson's correlation coefficient) is that it provides a quantitative evaluation of the evidence against the null hypothesis [47]. It also relies on a direct pairwise comparison between hypotheses based on a probability ratio, a quantity which is easy to interpret. Bayesian statistics are susceptible to the same caveats as the frequentist approach in cases of spatial autocorrelation: an illusion of more independent data points than there really are, which may for instance exaggerate the relationship between two variables if their spatial structures are aligned [48]. In a Bayesian analysis, different hypotheses are compared to identify which is most likely, and what values its parameters are most probable to take on. Bayesian linear regression is performed in two steps: model fitting and model comparison.
For a given hypothesis H, Bayes' theorem states: where p(H|D) = the probability H given the data (D), i.e., the posterior probability of the model; p(D|H) = the probability of the data (D) given H, i.e., the likelihood; p(H) = the probability of H, i.e., the prior probability of the model; p(D) = the probability of the data (D), i.e., the evidence.

Model Fitting: Bayesian Linear Regression
During model fitting, Bayes' theorem is used to infer the probability distribution of a given model's parameters based on the observed data (D). This is done for each hypothesis being considered. Equation (2) presents a linear model for a hypothesis H denoted as (α, σ), where the model parameters are coefficients α and standard deviation σ. It is given here for the ith sample of a dataset. The noise term η i in Equation (2) comes from a Gaussian probability distribution centered on zero with parameter σ as a standard deviation, i.e., η i ∼ N(0, σ), as shown in Equation (3). It is assumed to hold all of the measured data's variability.
where y i = the variable being modeled, x i = a vector representation of the model variables, α = a vector representation of the coefficients, η i = a noise term depending on model parameter σ, where σ = a model parameter.
Rearranging Equation (2) to express the noise term η i as the difference between observations y i and model results x i · α, we have: By combining Equations (4) and (3): Assuming measurements are independent, the probability p(D|H) of observing the data given model (α, σ) (i.e., the likelihood for a particular set of parameter values) can be developped as a product of the individual probabilities p(y i , x i |α, σ) for each of the N samples in the dataset: p(D|H), which is the quantity needed for model comparison, can now be computed by combining Equations (5) and (6): Equation (7) is computed over a range of discrete values for all model parameters. The resulting probability distribution includes p max (D|0, σ), the maximum probability in the distribution when all α coefficients are zero, i.e., the null hypothesis H 0 , and p max (D|α, σ), the maximum probability that can be found within the parameter space, i.e., the non trivial hypothesis H 1 . Other hypotheses can be constructed from the linear model by setting only some of the coefficients to zero. The parameter space should accomodate the full distribution of each parameter's marginal probability, i.e., p(α 1 ), p(α 2 ), etc. Examples for our model are shown in the supplementary ( Figure S7). The source code for our analysis is available on GitHub [49].

Model Comparison: Bayes Factor
Models are compared based on their probability given the data, i.e., p(H i |D). Although these probabilities cannot be computed directly, their pairwise ratio can be estimated based on each model's posterior probability p(D|H i ). In the absence of a priori information on which values the model parameters will take, we assume a uniform prior probability p(H) for every possible set of parameter values. The p(H i ) therefore cancel out, as well as the evidence p(D) which is the same regardless of the hypothesis. This pairwise ratio is called the Bayes factor K; it quantifies the evidence in favor of one of the hypotheses.
Considering two hypotheses, e.g., a non trivial hypothesis H 1 and the null H 0 and using Equation (1), the Bayes factor is therefore: If K = 1, the data offers no evidence against H 0 . A value of K = 10 0.5 is deemed sufficient to give substantial evidence against H 0 , while K = 10 and K > 100 respectively denote strong and decisive evidence [50]. We used the Bayes factor to compare a lineardependence hypothesis H i with a null hypothesis H 0 , where the parameter associated with the linear dependency is set to zero. We computed the Bayes factor using each hypothesis' maximum posterior probability:

Application to Snow Depth, Ice Thickness and SAR
We explored the dependence of the C-and X-band SAR HH log-scale backscattering coefficient on snow depth and ice thickness by assuming a linear relationship with either snow depth h s via parameter γ (H snow ), ice thickness h i via parameter δ (H ice ), or with both snow depth and ice thickness (H both ). Our approach to investigating the linear correlation between variables under a Bayesian framework is similar to the one proposed by Wetzels and Wagenmakers [47]. From previous studies, we expect γ < 0 and δ > 0. The duallinear-dependence hypothesis H both is therefore a function of an offset parameter σ HH 0 , linear terms γh s and δh i , and a noise parameter η σ , as shown in Equation (10). The null hypothesis H 0 is Equation (10) with both the γ and δ parameters set to zero; the snow-only hypothesis H snow is the same equation with only δ set to zero; the ice-only hypothesis H ice is Equation (10) (Figures S10 to S13).

Geary's C for Spatial Autocorrelation
We used Geary's C to investigate the presence of autocorrelation in our snow depth and ice thickness measurements. Real environments are structured by physical processes, such as currents and winds, which create gradients and patches [51]. In a variable such as snow depth or ice thickness, these features manifest as spatial autocorrelation -nearby measurements take values that are most or less similar than expected for a randomly assigned pair. This dependence violates the independence assumption central to many statistical tests such as Pearson's correlation coefficient. Autocorrelation will generally exaggerate the significance of such a test, for instance lowering its p-value [52]. Geary's C is a statistic developed to quantify spatial autocorrelation and for which a p-value may be evaluated for significance. Positive autocorrelation translates to a C value between 0 and 1, and negative autocorrelation produces values greater than 1; the no-correlation value is C = 1 [51].

Results
Here we present spatial distributions for snow depth and ice thickness measurements and an estimation of their spatial autocorrelation using Geary's C and seasonal median values for polarimetric parameters H-α and the HH backscattering coefficient, for all cases in the study. We then present the linear correlation coefficients between these variables and the results of the Bayesian linear regression analysis. Figure 3 highlights the spatial distribution of snow depth by showing end-of-winter standardized thickness measurements, color-coded according to their deviation from the mean, and the associated Geary's C value when its p-value was above 0.05. These results serve to evaluate the reliability of the statistical analyses presented in Section 3.3. In Salluit, snow accumulation is concentrated in the center of the study area in 2017, with a C value of 0.1 indicating strong positive autocorrelation, while being rather heterogeneous in 2018. In Deception Bay, the 2016 distribution is also heterogeneous. In 2017, snow accumulated to the north-west of the study area, around Moosehead Island, as reflected by a C value of 0.2. In 2018, accumulation appears to follow a gradient aligned to the north-east, with the deepest snow along the north-eastern shore and in front of the fjord's transverse arm. In Kangiqsujuaq, snow depth seems to increase along a south-east gradient for the three years, with higher accumulation in front of the community and along the south-eastern shore. When tested using the Shapiro-Wilk test, snow depth measurements rejected the normality hypothesis for three cases: Salluit in 2017 and Deception Bay in 2017 and 2018 (see Figure S3 in the supplementary materials). Figure 4 shows the spatial distribution of ice thickness for the same cases. Geary's C could not be reliably determined for any of the cases. Notable spatial structures in Salluit include thicker ice to the south-west of the study area both in 2017 and 2018. In Deception Bay, the ice was generally thicker along the south-western shore for all years, and thinner around Moosehead Island in 2017. In Kangiqsujuaq, features vary from year to year: thinner ice along the fjord's deepest area in 2016, a south-eastern gradient in 2017, and thicker ice in the north of the broader 2018 study area. When tested using the Shapiro-Wilk test, ice thickness measurements rejected the normality hypothesis for two cases, namely Deception Bay in 2017 and 2018 (see Figure S3 in the supplementary materials). Figure 5 shows seasonal median entropy (H) and alpha angle in the C-band, as well as HH backscattering coefficients in the C-band and the X-band when available, for each site and each season. All cases fall in the region of the H-alpha plane associated with surface scattering, i.e., alpha < 40 • [39]. For Salluit, 2015-2016 stands out with a higher alpha angle than the other seasons, suggesting a bigger volume contribution, and a lower C-band HH backscattering of −21 dB. Of all the sites and cases put together, Deception Bay 2015-2016 is the closest to the region of volume scattering. Its C-band backscattering is also the lowest with −25 dB compared to −17 dB in the next two years. In the X-band, the backscattering is also lowest that year with −20 dB. The difference in backscattering intensity between the C-and X-band is −5 dB for 2015-2016, zero in 2016-2017, and −1 dB in 2017-2018. For Kangiqsujuaq, the alpha angle was generally higher than for the other sites except in 2015-2016, where it was surpassed by the remarkably high alpha angle for Deception Bay.    [39]. Bottom: HH backscattering coefficients in the C-band (full markers) and X-band (empty markers).

Relationship between SAR and End-of-Winter Snow Depth or Ice Thickness
The linear correlation coefficient (r-squared) between the C-band HH backscattering coefficient and either snow depth or ice thickness at the end of the winter is between zero and 0.1 for all cases in Salluit. In Deception Bay, it ranged from zero to 0.3 for the C-band and zero to 0.6 for the X-band, and was always zero in Kangiqsujuaq ( Figures S4-S6 in the supplementary materials). We also investigated the following parameters: VV and HV backscattering coefficients, the VV/HH co-polarization ratio, the HV/HH and VH/VV cross-polarization ratios, and the entropy, anisotropy and alpha angle. They did not perform systematically better than the HH band as shown in the supplementary materials (Tables S3 and S4), despite the VV band sometimes having an r-squared 0.1 higher than the HH band or the alpha angle have a linear correlation coefficient of 0.5 with both snow depth and ice thickness in Deception Bay 2017. Figure 6 shows the most likely hypothesis from the Bayesian linear regression analysis of the relationship between the C-or X-band HH backscattering coefficient and either snow depth, ice thickness, or both. The associated pairise Bayes factor K is shown in the supplementary materials ( Figures S8 and S9). In the Salluit cases, the most likely hypotheses for 2017 and 2018 were the null and a negative relationship between C-band backscattering and snow depth, respectively. In Deception Bay, the most likely hypothesis for the C-band in 2016 was a negative relationship with ice thickness, while the null was most likely for the X-band. In 2017, it was a negative relationship with snow depth for both the C-and X-band. In 2018, the most likely relationship was a positive one with ice thickness for both bands. In Kangiqsujuaq, no hypothesis was more likely than the null for the three C-band cases.

Discussion
In this section, we discuss the spatial distribution of environmental variables ice thickness and snow depth, the relationship between C-or X-band backscattering and either variable, the interpretation of Bayesian linear regression results and perspectives for future work.

Spatial Structure in Environmental Variables
Here we investigate spatial structure in snow depth and ice thickness, in part because of its possible impact on our correlation analysis. Consider similar spatial structures in two environmental variables A and B (e.g., ice thickness and ice surface roughness), caused either by chance or by a physical process which affects both of them. An observed correlation between variable A and a third variable C (e.g., SAR backscattering) could be the accidental byproduct of a correlation between variables B and C, mirrored onto variable A because of its spatial structure's alignment with that of variable B.

Snow Depth Distribution
According to Geary's C, snow depth presented spatial autocorrelation in half of the cases, while the test was inconclusive for the other half ( Figure 3). A total of 30 measurement points was sometimes insufficient to get a conclusive statistic. This precluded the use of more advanced spatial structure analysis tools such as correlograms. We therefore conservatively assume all cases presented spatial autocorrelation.
Wind action is known to cause preferential snow accumulation patterns, such as in Van Mijenfjorden, a Norwegian fjord of the Svalbard archipelago [53]. The heterogeneous snow depth distribution observed in Salluit in 2018 (Figure 3) suggests along-fjord winds might have transported snow in and out of the study area. Spatial structure would then be caused by another factor than wind alone, such as ice surface roughness [6]. By contrast, in 2017 the snow accumulated in the middle of the Salluit study area. Because the C-band backscattering was similar in both seasons, suggesting comparable ice roughness, we speculate that dominant winds were oriented parallel to the fjord length during depositional or drifting events in 2017. There are no recent wind data available for Salluit; 2001 data [54] give northern and north-eastern winds as dominant. Since Inuit have reported more variable winds [2], readers should be aware that the situation might have changed in 20 years.
In Deception Bay, we suggest that heterogeneous snow depth in 2016 (Figure 3) might be associated with the remarkable smoothness of the ice that year, as illustrated by a C-band HH backscattering coefficient of −25 dB ( Figure 5) and documented in a previous study [31]. Snow accumulation around Moosehead Island in 2017 suggests a predominance of the along-fjord dominant winds, both north-western and south-eastern). In contrast, greater snow depth along the north-eastern shore and at the outset of the fjord's transverse valley in 2018 is rather consistent with the effect of across-fjord dominant winds. Indeed, transverse winds may be funnelled by valleys leading into the fjord [55]. We do not have up-to-date information on dominant winds in Deception Bay; the most recent observations date back to 1963-1973, more than 40 years ago [28]. In Kangiqsujuaq, snow depth consistently exhibited a south-eastern gradient (Figure 3). This is consistent with across-fjord dominant winds reported by NAV Canada in 2001 [54].

Ice Thickness Distribution
None of the cases included enough data to quantify spatial autocorrelation in ice thickness using Geary's C. Yet, visual interpretation of the ice thickness distributions (Figure 4) shows that the ice was generally thicker towards the south-western end of the study area in Salluit, where the water is more shallow. The ice was also thicker along the south-western shore in Deception Bay, where the snow is usually thin. In 2017 and 2018, ice thickness patterns in Deception Bay seemed to mirror those in snow depth, with thicker snow leading to thin ice and vice versa, similarly to observations made in Van Mijenfjorden (Svalbard, Norway) [53]. In Kangiqsujuaq, the ice in 2015 was thinner along the middle of the fjord's length, i.e., in the deepest part. It was thicker in the northern part of the broader study area in 2018, despite deep waters. We speculate that ice formation is more dynamic near the Hudson Strait and potentially included greater ice rafting and other accumulation processes. Near Pangnirtung in Nunavut, fjord outlets are known to present strong currents, even preventing solid ice formation [56].

Case A: Thin Snow Cover
Deception Bay 2018 is one of only two cases where the ice thickness hypothesis was most likely to explain C-or X-band HH backscattering, here through a positive relationship ( Figure 6). It also presented the strongest Bayes factor supporting a non trivial hypothesis and a linear correlation coefficient of 0.6 between ice thickness and the X-band (Figures S5, S8 and S9 and in the Supplementary Materials). This case is similar to other Salluit and Deception Bay cases in 2017 and 2018 in that it featured smaller alpha angles and higher backscattering ( Figure 5) than the smooth ice 2016 case in Deception Bay. Both indicators suggest surface scattering, which we attribute to slightly rougher ice due to freeze-up from nilas patches in dynamic conditions, as confirmed for one of the sites in a previous study [31]. It stands out from the other slightly rough ice cases with its thinner snow cover of eight cm on average (Table 3).
We speculate that in this case, (i) the brine-wetted snow was not thick enough to prevent the radar from reaching the ice surface, (ii) small-scale surface roughness allowed the surface scattering from the ice to make its way back to the radar, (iii) the dominant backscattering mechanism was surface scattering on saline ice and iv) the top-layer salinity, surface roughness or both were positively correlated with ice thickness which led to a positive correlation between the latter and HH backscattering. The Bayes factor associated with this correlation might also have been inflated by spatial autocorrelation, which we cannot rule out due to an inconclusive Geary's C (Figures 3 and 4). Note that the Deception Bay 2017 case also presented opposite snow depth and ice thickness gradients (Figures 3 and 4) but with thicker snow of 15 cm on average, and no correlation between backscattering and ice thickness was detected in that case. Our results for this thin snow case are similar to those of Nakamura et al. [24,25] in that we observed a positive correlation.

Case B: Thick Snow Cover
The three end-of-winter cases from Kangiqsujuaq stand out by their total lack of correlation between either snow depth or ice thickness and the HH backscattering coefficient ( Figure S6 in the supplementary materials). Similarly, none of the tested hypotheses was more likely than the null in the Bayesian linear regression analysis ( Figure 6). We attribute this behavior to the deeper snow found in Kangiqsujuaq compared to the other two sites: 22 to 25 cm on average during the study (Table 3). We speculate that in this case (i) the brine-wetted snow was too thick to allow the radar signal to reach the ice surface, (ii) the dominant backscattering mechanism was surface scattering on the interface between dry snow and brine-wetted snow.

Case C: Very Smooth Ice
Deception Bay 2016 is the second of only two cases where the ice thickness hypothesis was most likely in the Bayesian analysis, and only for the C-band ( Figure 6). The relationship is negative, contrary to case A for a thin snow cover. This case stands out due to a lower backscattering coefficient and a higher alpha angle than all others ( Figure 5). The ice cover was particularly smooth ice due to a thermal freeze-up [31] and its alpha angle ( Figure 5) indicates that volume scattering played a more prominent role that year for the C-band than in the other cases [39]. We made the same conclusion for the X-band data in a previous publication [31]. Our observed difference of 5 dB between the two bands ( Figure 5) is identical to reports by Nandan et al. in a similar case where volume scattering was deemed to be important [23].
In the C-band, the hypothesis for a negative linear relationship with ice thickness was most likely ( Figure 6). In the X-band, however, the null was most likely and the linear correlation coefficient between snow depth or ice thickness and either frequency's HH backscattering coefficient was zero ( Figure S4 in the supplementary materials). The relationship is therefore almost imperceptible and in the opposite direction than in case A for a thin snow cover over slightly rougher ice. We speculate that (i) the brine-wetted snow was not thick enough to prevent the radar from reaching the ice surface, (ii) small-scale surface roughness was too small for the surface scattering from the ice to make its way back to the radar. Note that what little relationship there is between ice thickness and backscattering is negative in this case, which differs from previous reports [24,25].

Relationship with Other SAR Parameters
The C-band co-polarization and cross-polarization ratios did not perform significantly better than the HH band. The strongest linear correlation coefficient we observed was 0.3 between the cross-polarization ratios and the ice thickness in the Salluit 2017 and 2018 cases ( Table S3 in the supplementary materials). This is different from results by Nakamura et al. at a comparable incidence angle of 37 • [57]. They reported an improved correlation with ice thicknesses up to 120 cm for the co-polarization ratio (r-squared = 0.6) compared to the correlation with the HH and VV backscattering coefficients themselves (r-squared = 0.4 and 0.3). The authors attribute this to ratio sensitivity to the dielectric constant of the toplayer ice [24], i.e., its brine content, and cancelling out of the small-scale surface roughness effect. The fact that the co-polarization ratio hides what little correlation there was between ice thickness and C-band SAR in Deception Bay 2018, with an r-squared of 0.0 compared to 0.3 (Table S3 in the supplementary), suggests that surface roughness was responsible for at least part of the observed correlation.

Snow Depth vs. C-or X-Band Backscattering
There are only two cases where the snow thickness hypothesis was most likely to explain the HH backscattering, both for the C-and X-band: Deception Bay in 2017 and Salluit in 2018 ( Figure 6). However, the associated linear correlation coefficients are very low, either 0.1 or 0.2 depending on the frequency ( Figures S4 and S5 in the supplementary materials). These cases are similar to a third, Salluit 2017, in terms of average snow thickness (Table 3), position in the H-alpha plane and median HH C-band backscattering ( Figure 5), and near-zero linear correlation coefficient between the two variables ( Figure S4 in the supplementary materials). We were therefore unable to reproduce observations by Gill et al. [27] of a moderate positive correlation (r-squared > 0.5) between snow depths ranging from 5 to 35 cm and C-band HH backscattering at comparable incidence angles of 32 • to 36 • .

Interpreting Results from a Bayesian Linear Regression
Despite very low linear correlation coefficients of zero to 0.3 between either thickness variable and the C-band HH backscattering coefficient ( Figures S4 to S6 in the supplementary materials), the Bayesian linear regression analysis gave a non trivial hypothesis as more likely than the null for half of the eight cases ( Figure 6). While these four cases all presented at least substantial evidence [50] against the null hypothesis (K > 10 0.5 ), the evidence in the Deception Bay 2018 case was decisive (K > 100) both with the Cand X-band ( Figures S8 and S9 in the supplementary materials). In that case, the linear correlation coefficients with ice thickness were 0.3 and 0.6 ( Figure S4 in the supplementary materials). This illustrates the need to contextualize results from a Bayesian model comparison analysis with a familiar indicator such as Pearson's linear correlation coefficient. The apparent disparity between the Bayesian and frequentist approaches, for example in case A, can be traced back to the fact that the Bayesian linear regression model evaluates the probability that a linear relationship exists considering the data provided. In contrast, the linear correlation coefficient provides a quantitative assessment of the degree to which the variables are correlated. The Bayesian hypothesis testing allowed us to identify the cases where the null hypothesis was more likely than a linear relationship with either thickness variable. For the other cases, we used it to identify which variable carried the most information about the backscattering coefficient.

Perspective for Future Work
Further work is needed to elucidate which of small-scale surface roughness or toplayer ice salinity is correlated with ice thickness in such a case. The freeze-up process and resulting small-scale roughness at each measurement location would be useful information for such a study, and processes such as the impact of snow redistribution on snow and ice salinity profiles should be taken into account. Surface desalination from snow redistribution might be different according to ice roughness and sheltering of the study area, and may differ between landfast and drift ice. Additionnaly, ice formed in an open area more than 250 km wide in the southern portion of the Sea of Okhotsk might be difficult to compare with ice formed in sheltered fjords less than 5 km wide, as evidenced by our lower X-and C-band HH backscattering coefficients. This might explain disparities between this study on fjord landfast ice and studies by Nakamura et al. on drift ice in the Sea of Okhotsk [24,25]. Future sampling should be designed to ensure that spatial autocorrelation may be quantified, for instance using Geary's C. Grid-like sampling may not be the best strategy to capture the spatial variability of snow depth and ice thickness. Local experts can help design a sampling strategy which would cover extremes and may help outline homogenous areas. Careful treatment of spatial autocorrelation is essential because it can inflate our confidence in the detection of a correlation between two variables [48,52]. While none of the polarimetric parameters we explored performed systematically better than the co-polarized backscattering coefficients, we did not examine the compact polarimetry parameter developed by Zhang et al. [26]. Following the launch of the RADARSAT Constellation Mission in 2019, access to RCM data started in 2020 and will continue to improve. Finally, although the incidence angles used in this study were already steep (35 • to 38 • ), it would be interesting to see if even higher incidence angles could improve our results. Indeed, Gill et al. [27] saw that the strength of the linear correlation between snow depth and C-band HH backscattering increased with incidence angle between 26 • -28 • to 35 • -37 • , and in the Nakamura et al. studies [24,25] reporting a strong correlation between ice thickness and HH backscattering, the incidence angles were 45 • and 39 • for the C-and X-band, respectively.

Conclusions
In this article, we combined field measurements and satellite SAR data in order to characterize seasonal snow-covered sea ice from 2015 to 2018 in three Nunavik fjords of the Hudson Strait.
In cases of landfast and undeformed first-year sea ice, we conclude that the necessary conditions for detecting a correlation between sea ice thicknesses above 30 cm and SAR HH backscattering in the C-or X-band are a snow cover thinner than 10 cm and slightly rough ice formed from dynamic processes. In terms of usability for predictions, the correlations we observed in these conditions were poor at best for the C-or X-band, respectively, with r-squared values of 0.2 and 0.6. In cases with snow thicker than 20 cm on average or with very smooth ice, no correlation could be detected with ice thickness. No correlation above 0.3 was observed between snow depth and backscattering. The Bayesian linear regression analysis proved to be useful in categorizing each case according to their most likely hypothesis out of the ones we tested. Our results differ from previous reports of a correlation either between ice thickness and the co-polarization ratio in the C-band [25] or between snow depth and the HH backscattering coefficient in either frequency [27].
Backscattering in the C-and X-bands was either different or similar depending on the type of ice. Over smooth ice formed from thermal freeze-up, backscattering from both bands has a significant volume scattering contribution and their HH backscattering coefficients present a 5 dB difference. Over slightly rougher ice formed from consolidated nilas patches, surface scattering dominates at both frequencies. Their HH backscattering coefficients differ by 1 dB or less.
Supplementary Materials: The following are available at https://www.mdpi.com/2072-4292/13/4 /768/s1, Figure S1: Monthly total precipitation measured at Salluit airport, Figure S2: Bathymetry for Salluit, Deception Bay, and Kangiqsujuaq, Figure S3: Measurement distributions for snow depth and ice thickness measured in April-May of each season, Figure S4: C-band backscattering HH coefficient versus snow depth and ice thickness for available cases in Salluit, Figure S5: Same as S4 for Deception Bay, including X-band data, Figure S6: Same as S4 for Kangiqsujuaq, Figure S7: Marginal probability distributions for the H snow hypothesis parameters in the Salluit 2018 case, Figure S8: Pair-wise logarithmic Bayes factors for the relationship between the C-band HH backscattering coefficient and either snow depth, ice thickness, or both, Figure S9: Same as S8 for the X-band in Deception Bay, Figure S10: Mean and standard deviation of the parameters' marginal probability distributions for H snow applied to C-band data, Figure S11: Same as S10 for H ice , Figure S12: Same as S10 for H both , Figure S13: Same as S10 for H snow , H ice , and H both , applied to X-band data; Table S1: Fieldwork dates in Salluit and Kangiqsujuaq, RADARSAT-2 acquisition dates, and the difference between the two, Table S2: Fieldwork dates in Deception Bay, RADARSAT-2 and TerraSAR-X acquisition dates, and the difference between the two, Table S3: Linear correlation coefficient (r-squared) between C-band SAR parameter values and April-May snow depth and ice thickness.

Data Availability Statement:
The data presented in this study will be available in the public repository PANGAEA in 2021. They include : (i) SAR parameter statistics at AOIs for every RADARSAT-2 and TerraSAR-X image, (ii) SAR parameter values at snow or ice thickness measurement locations for each case, and (iii) shapefiles of snow depth and ice thickness measurements for each case [58]. The code used to process the RADARSAT-2 images using ESA's SNAP is available at https://gitlab. com/sdufourbeausejour/java-snap, accessed on 12 February 2021 [35]. The code used to compute SAR parameter median from AOIs is available at https://github.com/sdufourbeausejour/tiffstats, accessed on 12 February 2021 [45]. The code used to extract SAR parameter values at shapefile feature locations is available at https://github.com/sdufourbeausejour/tiff_at_shp, accessed on 12 February 2021 [46].