SAR data and processing

Details of the SAR scenes used in this study are provided in Table 1. It is important to note that the portions of the Sentinel-1 SAR scenes that are spatially coincident to high-resolution VIS-NIR scenes from which f _p data cover a narrow range of radar incidence angles (θ). The archived ASAR data from AB were selected to enable us to further investigate the role of θ on backscatter and f _p relationships.

Table 1. Basic image characteristics of the SAR scenes acquired for this study

Pass refers to ascending (Asc) and descending (Des) orbits; Swath to the available swaths of Envisat-ASAR APP mode products (IS1–IS7) and Sentinel-1 products; and incidence angle (θ) ranges of the overall swath and for the portion examined.

Three Sentinel-1 Extra Wide Swath Mode (EW) scenes, S1–S3, were acquired for sites VN and VS. The EW mode was designed for maritime use, and particularly for imaging sea ice, by employing a wide swath. In EW mode, Sentinel-1 acquires data over a 400 km swath at 20 m by 40 m spatial resolution (ESA, 2013). Data are acquired in one of four independent polarization configurations using one H or V polarization transmit chain, and one or two parallel receive chains for H and V polarizations. This produces either a single polarization VV or HH, or a dual-polarization HH + HV or VH + VV, scene comprised of five sub-swaths and spanning a θ range of 19–47°. Scenes S1–S3 were processed to medium resolution level-1 ground range detected format with 40 m by 40 m pixel spacing, prior to downloading from the Data Hub (ESA, 2013). Sentinel-1 is an ongoing constellation mission with two satellites, Sentinel-1A and Sentinel-1B, sharing the same near-polar orbital plane with a 180° orbital phasing difference (ESA, 2013). Sentinel-1A was launched in early 2014 and Sentinel-1B in early 2016.

Three ASAR Alternating Polarization Mode Precision image (APP) products of site AB, named AB1–AB3, were acquired from the ESA Earth Observation Link. In APP mode, ASAR interleaved looks along-track in H and V polarizations within the synthetic aperture, resulting in the production of two ground range detected image channels (i.e. dual polarization images) out of three possible transmit–receive polarization combinations (HH + VV or HH + HV or VV + VH) at a nominal resolution of 30 m. The sensor was capable of acquiring one of seven different image swaths spanning a total θ range of 15–45°. Swath widths ranged from 56.5 to 104.8 km, depending on which of the seven swaths was chosen. The Envisat satellite was launched in early 2002 and ceased operating in early 2012.

SAR scenes were processed according to the flowchart in Figure 2. The thermal noise removal was done by subtracting the noise estimate values, retrieved from the annotation dataset of each SAR image product, from the corresponding image channel. Image bands were calibrated to the common sigma-nought backscatter σ _hh^° or σ _hv^° and speckle filtered using the Lee Filter and a 5 by 5 sliding window. A cross-polarization ratio R _hv/hh (σ _hh^°/σ _hv^°) band was then derived. Image bands were map projected to the WGS 1984/UTM projection using the Bilinear resampling method and clipped to the appropriate site extent, as determined by the coverage of the corresponding high-resolution VIS-NIR scenes at sites AB, VN and VS. Envisat-ASAR scenes were projected with 12.5 m pixel spacing, whereas Sentinel-1 scenes were project with 40 m pixel spacing. After clipping, an 8-bit scaled version of the σ _hh^° band from each SAR scene was exported for segmenting each site into discrete image objects for the object-based image analysis (described below).

Fig. 2. Flow chart depicting the SAR image processing chain. The final GeoTIFF image product contains calibrated backscatter sigma-nought (σ°), backscatter ratio and texture bands.

Second-order texture parameters were derived from each of the three clipped image bands from the six SAR scenes, σ _hh^°, σ _hv^°, and R _hv/hh, using the GLCM method developed by Haralick and others (Reference Haralick, Shanmugan and Dinstein1973). The co-occurrence parameters are based on computing the joint probabilities of all pair-wise combinations of brightness values in a spatial window (kernel) according to two parameters: a displacement (inter-pixel distance) and, an orientation (0, 45, 90 and 135°) (Clausi, Reference Clausi2002). From Clausi (Reference Clausi2002), this is defined as:

(1)$$\Pr (x) = \left. {\left\{ {C_{ij}} \right.\left \vert {\left( {\delta, \theta _{\rm o}} \right)} \right.} \right\},$$

where C _ij is the co-occurrence probability between grey levels i and j, δ is the inter-pixel distance, and θ _o the orientation. C _ij is defined as:

(2)$$C_{ij} = \displaystyle{{P_{ij}} \over {\sum\limits_{i,j = 1}^G {P_{ij}}}}, $$

where P _ij represents the number of occurrences of grey levels i and j within the kernel, given a δ and θ _o pair, and G is the quantization level. The texture parameters are then derived from statistics applied to the co-occurrence probabilities. The comprehensive review of the use of co-occurrence texture parameters for the separability of sea-ice types and sea-ice classification by Clausi (Reference Clausi2002) was used to guide the selection and calculation of texture parameters. Co-occurrence probabilities were computed for a 5 by 5 sliding window, using G = 64, and values derived from all θ _o (0, 45, 90 and 135°) were averaged. For δ, a value of 2 was chosen to negate the effect of correlated neighboring pixel values on texture calculations. As with SAR images in general, the scenes used here were oversampled, so that pixel spacing is approximately one-half the resolution in range and azimuth (ESA, 2013). SAR texture features contrast (CON), homogeneity (HOM), energy (ENE), entropy (ENT) and GLCM variance (GLV) were calculated (Table 2).

Table 2. Co-occurrence based texture statistics

VIS-NIR data and processing

A Quickbird image product with four channels spanning 450–900 nm and 2.4 m ground sample distance (GSD) and a panchromatic channel with 0.61 m GSD covered a 64 km² area of smooth FYI at site AB on 26 June 2006 at 18:44 UTC. Approximately 15% of the scene area that was cloud and land contaminated was manually masked out. Two GeoEye-1 image products containing four channels spanning the VIS-NIR range (450–920 nm) with a 1.7 m GSD, and a panchromatic channel with 0.41 m GSD covered 72 and 118 km² areas of mainly FYI, and mainly MYI, respectively, at site VN on 25 June 2015 at 18:20 UTC. A GeoEye-1 image product with the same bands covered a 100 km² area containing FYI and DFYI at site VS on 21 June 2016 at 18:45 UTC.

The high-resolution images were map projected to the WGS 1984/UTM projection. The Gram Schmidt pan-sharpening algorithm was used to fuse VIS-NIR bands from each product with its panchromatic pair, producing an output three-band (RGB) pan-sharpened image with a cell size of 0.5 m. The supervised maximum likelihood approach was used to partition the pan-sharpened scenes at AB and VN into binary classified images composed of snow/ice (0) and melt pond (1). A third class, drained melt pond (3) was added to the classification of the pan-sharpened scene at VS. By adding this class, we found that we could improve overall classification accuracy while still merging the two melt pond-related classes into one f _p class that is more representative of conditions prior to pond drainage. We found this was necessary since the VS scene was captured almost 20 d after pond formation, so drainage from FYI was evident in the scene. Training data for the supervised classification of each image were derived from 100 randomly chosen pixel locations. An accuracy assessment was also performed on each image, using 100 randomly chosen pixel locations for comparisons between the classification output and the input image. The spectral differences between light snow/ice patches and dark melt ponds made for relatively straightforward classifications, and overall accuracies were >97%. Examples are shown in Figure 3.

Fig. 3. Subsets of 300 by 300 m size showing true-color representations of surface conditions corresponding to the three main ice classes investigated in this study (top) and classification results with ice colored white, melt pond colored dark blue and drained melt pond colored light blue (bottom). (a) FYI in VN; (b) MYI in VN; (c) DFYI in VS.

Object-based image analysis

One drawback of using EO data is that it contains an arbitrarily defined, uniform sampling grid which is usually not appropriate for representing the entities of interest. In addition, the chosen scale for representing phenomena will affect results and interpretations derived from analyses. This is commonly referred to as the scale problem (Marceau, Reference Marceau1999). An object-based image analysis aims to mitigate the scale problem by segmenting the EO image into objects or segments that represent real-world entities in a manner consistent with patches in the landscape ecology literature (Turner and others, Reference Turner, Gardner and O'Neill2001; Benz and others, Reference Benz, Hofmann, Willhauck, Lingenfelder and Heynen2004). Here the regions of interest are constrained by the extents of the high-resolution VIS-NIR image products obtained for sites AB, VN and VS. The objects within each region of interest are individual pans or floes of sea ice which have unique dynamic and thermodynamic histories and, in a spatial context, have internally coherent and externally heterogeneous geophysical properties such as roughness which means they are: (a) likely to be unique areas in terms of thermodynamic evolution, hence winter state to spring melt pond formation; and (b) distinguishable in winter C-band SAR images.

The 8-bit scaled σ _hh^° image bands produced during SAR processing (see Fig. 2 and Table 1) were segmented into image objects using the bottom–up, region-merging technique detailed in Benz and others (Reference Benz, Hofmann, Willhauck, Lingenfelder and Heynen2004) and implemented in the eCognition^® software. This algorithm was chosen since it allows the user to specify weightings that are used as homogeneity criteria for growing and segmenting coherent objects during the segmentation process. A color versus shape criterion, where color = 1 − shape, enables weighting toward color (i.e. backscatter intensity scaled to 8-bit) versus shape (i.e. spatial) information. The algorithm also allows the user to vary the sizes, and therefore the number of objects, created during the segmentation by specifying a scale parameter. User specified intensity and shape criteria are input into a heuristic, region growing algorithm which grows regions until the criteria are met. This results in a set of image objects of heterogeneous intensity and spatial characteristics. Additionally, a scale parameter can be used to create a multiresolution segmentation with multiple layers of image objects with successively larger (or smaller) scale.

A hierarchical segmentation approach was used to address the fact that sea-ice zones are not necessarily distinctly bounded, and no segmentation represents a perfect representation of distinct real-world entities. Three levels of segmentation were created by varying the scale parameter of the segmentation algorithm. The hierarchical approach also enabled us to assess the impact of object scale and resultant data aggregations on the strengths of correlation coefficients and coefficients of determination. The three segmentation levels created image object sets composed of progressively larger and fewer image objects for overlaying on spatially coincident SAR and optical scenes. We named the object sets fine, sharp and coarse and calculated the areas of the objects making up each of those sets (Fig. 4).

Fig. 4. Image object sets for zone of predominantly multiyear sea ice in Victoria Strait north (VN) collected in 2015. Relatively smooth first-year sea ice is evident in the lower portion of the scene. (a) Eight-bit scaled HH polarization band of Sentinel-1 image VN1 (see Table 1) used for segmentation into three image object sets comprising progressively smaller and fewer image objects from (b)–(d). The coarse segmentation produced n = 42 objects (b). The sharp segmentation produced n = 87 objects (c). The fine segmentation produced n = 220 objects (d).

From fine to coarse scales, the object sets ranged from n _obj = 237 to n _obj = 37 (e.g. see Fig. 4). Average object sizes at the fine scale ranged according to the object set (AB, VN, VS) from 0.21 to 0.57 km². At the sharp scale, the range was 0.56–1.1 km², and at the coarse scale 1.3–2.8 km². Each object sampled on average 210–355 Sentinel-1 pixels at the fine scale, 493–817 at the sharp scale, and 1057–1392 at the coarse scale (VN and VS). The smaller 12.5 m pixel spacing of the ASAR images meant that on average 1344–8685 pixels were sampled by objects at fine to coarse scales.

Data analysis

Objects from each site were labeled according to the dominant ice type contained within their bounds by using the winter SAR and spring VIS-NIR imagery, as well as ancillary information such as ice charts, to guide the analysis. Objects from AB all contained FYI. Objects from VN were labeled either FYI or MYI, since the two VN scenes comprised either smooth, thermodynamically grown FYI or MYI floes. The MYI-dominant scene from VN is shown in Figure 5, which also illustrates the segmentation and labeling process. Objects from VS were labeled FYI or DFYI, with the latter evident as predominantly linear and brightly backscattering features in SAR imagery (Fig. 6).

Fig. 5. (a) A 20 by 20 km subset of a Sentinel-1 SAR scene of VN containing FYI and MYI; (b) The same scene with the object set overlaid on it; (c) The spring VIS-NIR scene of the same area; and (d) Result of the labelling of objects as FYI (black) or MYI (white).

Fig. 6. (a) A 20 by 20 km subset of a Sentinel-1 SAR scene of VS containing FYI and DFYI; (b) The same scene with the object set overlaid on it; (c) The spring VIS-NIR scene of the same area; and (d) Result of the labelling of objects as FYI (black) or DFYI (white).

Each of the three object sets fine, sharp and coarse was overlaid on the SAR scenes of calibrated backscatter and texture. Object metrics were computed from individual objects as the mean of pixel values within an object. Each of the three objects sets pertaining to a VIS-NIR derived snow/ice and pond binary classification output was used to calculate the f _p for each image object. This enabled an output database of winter backscatter and texture parameters, and spring f _p, for each of the scenes in Table 1 for statistical comparisons.

An exploratory data analysis was used to characterize the backscatter and f _p behavior of the sea-ice cover, and to assess backscatter levels relative to SAR sensor noise levels. A correlation analysis of object metrics enabled the examination of inter-relations between SAR backscatter, SAR texture and f _p. The correlation analysis focused on the linear association between backscatter and f _p values at fine, sharp and coarse scales. Assessment of the three ASAR scenes AB1–AB3 allowed us to examine the combined roles of object size (scale of aggregation) and θ on backscatter relationships with f _p. Linear associations between texture and f _p were calculated for the intermediate scale sharp image object set only. The Pearson's product–moment correlation coefficient (r) was used as a measure of linear association. The conventional deciBel (dB) transformation was applied to backscatter data (10*log10(σ°)) and, in some instances, the natural log transformation was applied to texture (log₁₀) to achieve linear association between variables and f _p.

Linear regression models for predicting f _p from winter Sentinel-1 backscatter were constructed using data from the intermediate scale sharp image object set. Only backscatter and GLCM texture parameters derived from the HH polarization channel were used, and input data from only one site were used in order to control for temporal fluctuations in pond fraction (see Melt pond prediction model section). All texture parameters were used as potential explanatory variables, with log-transformed texture parameters used where appropriate.

The backward elimination method was used to develop regression models. Backward elimination is a stepwise method where all possible predictor variables that pass a multicollinearity tolerance criterion are entered into the model, and highly multicollinear variables are removed. Following entry, independent variables with the smallest partial correlation with the dependent variable are sequentially passed through stepping method criteria. If the variable does not contribute a significant increase in r ² (α > 0.1) it is removed. This process is continued until the model with the smallest number of independent variables remains.

Site conditions

Field data were collected as part of a dedicated sea-ice geophysical and remote-sensing study at the AB site from 5 May to 15 July 2006. In 2016, a multidisciplinary sea-ice campaign was conducted on smooth FYI at DS near Cambridge Bay from 15 April to 22 June, though this location was ~150 km west of site VS. We chose site VS for some remote-sensing studies during the 2016 campaign since it provided a diversity of ice conditions lacking at the field site. No field data were collected in 2015.

Air temperature data were used to guide the selection of ASAR and Sentinel-1 scenes originally acquired during winter (pre-melt) conditions. In the seasonal period leading up to the ASAR acquisitions at AB, the maximum hourly 10 m air temperature was −9.6°C. The ASAR images of AB were collected at hourly air temperatures of −14.6°C or lower. The maximum hourly 10 m air temperature at the Environment Canada Cambridge Bay airport station (69.11N, −105.14W) in the period leading up to the Sentinel-1 acquisitions of VN in 2015 was −13.2°C. The Senitnel-1 images of VN were collected at an hourly temperature of −26.9°C. In 2016, air temperatures recorded at the same station showed a maximum temperature of −15.2°C in the lead up period and −27.2°C at acquisition.

Figure 7 shows time series f _p measurements made at site AB, within the same area covered by the acquired VIS-NIR scene from that campaign, and at site DS, 150 km east of the VIS-NIR scene acquired at VS during that campaign. Field measurements of f _p came from manual transects, where a surface class and depth (not shown) was recorded every 0.5 m along a 200 m line perpendicular to the major axes of melt ponds and snow/ice patches. Recorded classes were: snow, bare ice, melt pond, frozen melt pond/ice lid and slush/mixture; 400 points were logged during each transect. This approach is reliable for recording an areally representative f _p on smooth FYI only, i.e. when there is a quasi-uniform distribution of ponds and ice patches determined by the pre-melt alignment of snow dunes (Grenfell and Perovich, Reference Grenfell and Perovich2004; Istomina and others, Reference Istomina2015). It is a useful approach for recording the time series evolution of f _p and other surface features, though consideration must be given to estimation error due to mis-alignment. We tested for mis-alignment error by plotting the location of the in situ transect line on the classified VIS-NIR scene from AB, and determining the f _p for that line and lines at 10° intervals to ±20°. The f _p range across all angles was 0.06. Figure 7 also shows the timings of VIS-NIR acquisitions at AB and VS, along with means and SDs of f _p derived from classified scenes.

Fig. 7. Time series melt pond fraction recorded at the Allen Bay (AB) site in 2006, and at the Dease Strait (DS) site in 2016, from the date of pond onset. Orange and blue markers denote the acquisition times of high-resolution VIS-NIR scenes capture at AB in 2006, and Victoria South (150 km east of DS), respectively. The markers also indicate the means and SDs of melt pond fraction derived from the VIS-NIR scenes after classification.

From Figure 7, it is apparent that the VIS-NIR scene at AB was acquired 3.8 d after melt ponds formed, when f _p was rapidly increasing but still <0.5. This timing is just prior to the seasonal f _p peak for smooth FYI, when f _p is growing to its seasonal maximum due to a rapid flux of melt water from the decaying meteoric snow cover (Eicken and others, Reference Eicken, Krouse, Kadko and Perovich2002). When the VIS-NIR scene at VS was acquired, our observations on smooth FYI at DS indicated that ponds formed 19.8 d earlier. The f _p at DS was high at ≥0.7, though the snow cover had ablated and fluctuations in f _p linked to diurnal melt rates were evident. There is a good agreement between field and satellite-derived f _p at AB, whereas an inconsistency between field-derived f _p from DS, and satellite-derived f _p from VS, is apparent. Note that the offset in f _p between DS and VS is due to the inclusion of DFYI and FYI at VS, whereas at DS only FYI is captured (see Fig. 6). Since we did not have field data linked to the VN acquisitions, an examination of cloud-free MODIS reflectance images indicated that the VIS-NIR images were acquired 14 d after pond onset.

Figure 8 shows f _p distributions according to major ice types FYI and MYI derived from all acquired VIS-NIR scenes. From Figure 8, it is evident that a wide range of f _p on FYI and on MYI is presented by the VIS-NIR image data. As such, a wide range of ice conditions are included in the comparison to backscatter and texture statistics. The average f _p for FYI image objects ranged from 0.07 at AB to 0.91 at VS. The lowest mean f _p of 0.39 recorded at AB was measured shortly after pond formation began and f _p was still increasing. The FYI f _p distributions for VN and VS are very similar and representative of conditions after melt ponds have fully developed. Overall, the f _p statistics are similar to other observations of FYI and MYI made at a similar scale, e.g. from aircraft observations (Perovich and others, Reference Perovich, Tucker and Ligett2002; Eicken and others, Reference Eicken, Grenfell, Perovich, Richter-Menge and Frey2004).

Fig. 8. Melt pond fraction statistics derived from classified high-resolution VIS-NIR images, using the sharp image object set at each site to calculate melt pond fractions.

Relationship between winter backscatter and spring melt pond fraction

Incidence angle and scale: Allen Bay (AB) site

Correlations between backscatter parameters σ _hh^°, σ _hv^°, and R _hv/hh, and spring f _p for all three ASAR scenes at AB and all three aggregation scales, are given in Table 3. This enables a detailed examination of the combined roles of θ and scale (object size) on backscatter-f _p relationships. The correlations in Table 3 are all significant at α = 0.01 (P-value <0.01), apart from R _hv/hh and f _p from AB2 derived using the coarse segmentation (r = 0.362). The polarization ratio R _hv/hh and f _p from scene AB1 are strongly positively correlated at coarse (r = 0.846) and sharp (r = 0.838) scales. This relationship points to a strong sensitivity of R _hv/hh to the winter ice structure that leads to f _p formation on the relatively smooth FYI forming in the sheltered bay at AB. This relationship is strongest for AB1, where backscatter was measured at shallow θ, here 39.1–42°. Relationships between R _hv/hh and f _p are much weaker at steeper θ.

Table 3. Correlations r between Envisat-ASAR derived σ _hh^°, σ _hv^° and R_hv/hh, and melt pond fraction for FYI samples taken from all three AB scenes AB1–AB3, and by using all three scales of aggregation coarse, sharp and fine

In Table 3 σ _hh^° and σ _hv^° are negatively associated with f _p and the strongest correlations are found using the sharp object set. Generally strong correlations between σ _hh^° and f _p are found in data obtained from all scenes spanning the near to far θ range examined. Focusing on results in Table 3 derived using the sharp object set only, the strongest relationship between σ _hh^° and f _p comes from the shallow θ. The strongest relationship between σ _hh^° and f _p is found from scene AB2 (r = −0.792), acquired at the steepest θ of the three AB scenes. As near-range σ _hv^° is more strongly correlated with f _p than it is at mid- to far ranges, this implies that near-range σ _hv^° is more closely related to variations in surface structure on relatively smooth FYI that also lead to variations in f _p shortly after pond onset. Given the low levels of σ _hv^° intensity from smooth FYI, caution is needed when interpreting σ _hv^° or R _hv/hh. We found the overall σ _hv^° range from the three ASAR scenes to be −18.9 to −23.6 dB. Scenes AB1 and AB2 (ASAR swaths IS6 and IS2, respectively) had σ _hv^° consistently above the noise-equivalent sigma-zero (NESZ) of those swaths. Scene AB3 (ASAR swath IS4) σ _hv^° measurements were below the NESZ of that swath (ESA, 2006). The relatively consistent behavior of the σ _hh^° across the θ range investigated, combined with high σ _hh^° intensity levels relative to the system NESZ (here −9.2 to −20.4 dB) suggests the parameter can be used for the development of f _p retrieval techniques spanning a large portion of the wide 400 km swath employed by Sentinel-1 operating in EW mode.

Correlations between GLCM texture parameters and spring f _p for all three AB scenes, aggregated using the sharp scale, is shown in Table 4. All variables in Table 4 are significantly correlated with f _p at α = 0.01 (P-value <0.01), though the log transformation is required in the case of CON. As before it must be stated that GLCM parameters derived from σ _hv^° or R _hv/hh from AB3 may be contaminated by system noise. The potential utility of GLCM parameters derived from the σ _hh^° to be used along with σ _hh^° in the development of a f _p estimation approach is apparent. All GLCM parameters derived from σ _hh^° at all examined θ, i.e. from AB1 to AB3, are strongly correlated with f _p provided the log transformation is used where appropriate (see Table 4). GLCM parameters from σ _hh^° are generally most strongly associated with spring f _p when measured at the shallow θ, here 39.1–42°.

Table 4. Correlation r between Envisat-ASAR derived GLCM texture parameters contrast (CON), homogeneity (HOM), energy (ENE), entropy (ENT) and GLCM variance (GLV), and melt pond fraction for FYI samples taken from all three AB scenes AB1–AB3

*Significant at the α = 0.05 level.

The prefix log is used to denote instances where the texture parameters were log scaled after they were derived from HH, HV and HV/HH bands. All correlations are significant at α = 0.01 unless otherwise noted. Bolded text is used to highlight situations where the log transformation of a parameter improved the strength of correlation with f _p.

Sentinel-1 backscatter and FYI, DFYI and MYI pond fractions

Correlations between winter Sentinel-1 backscatter parameters σ _hh^°, σ _hv^° and R _hv/hh, and spring f _p for FYI, DFYI and MYI are shown in Figure 9. Data in Figure 9 were produced from datasets collected at VN and VS, with backscatter and f _p aggregated at the intermediate (sharp) scale. Table 5 provides additional information on significance testing as well as correlation results obtained by using the two other scales for aggregating data. The two-tailed significance of each correlation was tested at α = 0.01 and all correlations are significant (P-value = 0.00). The scatterplots in Figure 9 aid the interpretation of the winter backscatter behavior as it relates to spring f _p, as well as intensity levels relative to the Sentinel-1 NESZ. FYI and MYI are compared separately from FYI and DFYI, given that those datasets were collected at different locations and in different conditions.

Fig. 9. Correlations between pre-melt σ _hh^° and spring melt pond fraction (top), σ _hv^° and spring melt pond fraction (middle) and the cross-polarization ratio R _hv/hh and spring melt pond fraction (bottom) for FYI and MYI sampled at site VN (left panel). Correlations between pre-melt σ _hh^° and spring melt pond fraction (top), σ _hv^° and spring melt pond fraction (middle), and the cross-polarization ratio R _hv/hh and spring melt pond fraction (bottom) for FYI and DFYI sampled at site VS (right panel). All backscatter parameters are expressed in deciBel (dB) format. The Pearson's product moment correlation coefficient r is shown in each plot.

Table 5. Correlations r between winter Sentinel-1 derived σ _hh^°, σ _hv^° and R_hv/hh, and melt pond fraction for all FYI and MYI samples from site VN, and all FYI and DFYI samples from site VS. Testing was done using object sets coarse, sharp and fine in each case. P-values from two-tailed significance tests (t-tests) of each correlation are all 0. The number of sampled objects (n) from each object set are also given

FYI and MYI are easily distinguishable during winter due to the higher bulk salinity of FYI compared with MYI, which acts to prevent the penetration of C-band microwaves into its volume. Typical ice salinities for FYI are 5–8 ppt at 1–2 m thickness and 0.1–3 ppt for MYI (Weeks, Reference Weeks1981). MYI has much lower salinity due to desalination from a previous summer, thereby promoting the penetration of C-band microwaves and significant volume scattering from air bubbles and other constituents in the desalinated upper ice layer (Hallikainen and Winebrenner, Reference Hallikainen, Winebrenner and Carsey1992). This volume scattering mechanism enhances the backscatter of energy to the SAR and increases intensity levels relative to FYI. As well, thermodynamically grown FYI is relatively smooth compared with MYI which has meter-scale variations in topography due to the ablation and refreezing of melt ponds and hummocks from a previous year (or years’) melt. The large-scale variations in topography also enhance backscatter by increasing the occurrence of features titled toward the SAR, which enhance microwave backscatter. The σ _hh^° and σ _hv^° plots that combine FYI and MYI in Figure 9 show a continuum of backscatter intensity related to those ice types that is also connected to the ice topography, which controls the formation of spring melt ponds.

Comparing FYI and DFYI in Figure 9, there also appears to be a continuum of σ _hh^° and σ _hv^° intensity levels which are also connected to ice topography and spring f _p. As winter σ _hh^° and σ _hv^° increases due to enhanced backscatter from the ridges and ice blocks that make up DFYI, the spring f _p decreases due to topographical controls on melt water distribution from the same ridges and ice blocks (Eicken and others, Reference Eicken, Krouse, Kadko and Perovich2002). The σ _hh^° channel is more sensitive to this relationship; the σ _hv^° to f _p relationship is much weaker.

As with the Envisat-ASAR examined above, consideration must be given to low σ _hv^° measured by Sentinel-1 relative to the NESZ for the sensor and image mode. The ground range detected EW product used here has a maximum NESZ of −22 dB, though the σ _hv^° data in Figure 9 suggest sensitivity to sea-ice backscatter mechanisms down to −30 dB or less. Nonetheless, uncertainty regarding the utility of σ _hv^°, and by extension R _hv/hh, is exemplified by looking at the general behavior of those parameters in Figure 9. In the right panel of Figure 9, FYI and DFYI σ _hh^° decreases with increasing f _p, whereas σ _hv^° is almost static at ~−30 dB. The association between R _hv/hh and f _p switches from negative to positive in the left and right panels of Figure 9, despite the consistent behavior of σ _hh^°. Given that σ _hv^° levels are so low, combined with the fact that the HV polarization channel of Sentinel-1 is subject to variable swath-dependent noise levels across the sub-swaths which make up the EW and IW modes, backscatter and texture information from the HH polarization channel from Sentinel-1 should only be used.

Table 5 reveals generalizations to larger aggregation units that are matched by increases in the magnitudes of r. This is expected based on common geographic theory concerning modifiable areal units, where an increase in the magnitude of r between two variables occurs when those two variables are aggregated by larger units such as grid cells (Fotheringham and Wong, Reference Fotheringham and Wong1991). However, aggregating by larger units does not reveal the concomitant loss of variability in the dataset. Should a regression model be developed from the data created at the coarse resolution, with fewest aggregation units (n = 37 or n = 91), the explanatory power of the model would be expected to be weaker than would a model derived from a greater number of smaller aggregation units, in this case image segments.

Scatterplots showing relationships between winter GLCM texture parameters derived from σ _hh^° and spring f _p are given in Figure 10. Data were aggregated at the sharp scale, and log transformed texture parameters were shown in cases where the strength of association was stronger. All GLCM parameters in Figure 10 are significantly correlated with f _p α = 0.01 (P-value = 0.00). The texture variable logENE is not shown in Figure 10 since it is co-linear with HOM. Relationships between winter GLCM texture parameters and spring f _p are intuitive when considered in the context of the heterogeneity of spatial patterns of single-polarization backscatter intensity in SAR images as they relate to ice type. Relatively smooth FYI, which evolves to have high spring f _p, is also characterized by spatially homogeneous zones of backscatter intensity. MYI, which has a more complex topography and evolves to have lower f _p compared with FYI, is characterized by heterogeneous distributions of scattering mechanisms and pixel intensity values (see Fig. 5a). DFYI, which can have similar σ _hh^° intensity levels as MYI, also have more complex textures when compared with FYI. As such the texture parameters which are responsive to spatial heterogeneity or disorder, being logCON, ENT and GLV, are highest for winter MYI and DFYI and evolve to lower f _p. Texture parameters which are responsive to orderliness or uniformity, such as HOM and logENE (not shown), are highest for winter FYI evolve to higher f _p. This behavior of GLCM texture parameters of winter sea ice is not unfounded, as the inclusion of GLCM texture parameters has been shown to improve the classifications of sea-ice types in winter SAR images (Clausi and Deng, Reference Clausi and Deng2003). However, the further connection of textural features of ice types to spring f _p has not yet been made.

Fig. 10. Scatterplots of winter GLCM texture parameters and spring melt pond fraction for sites VN (FYI and MYI) and VS (FYI and DFYI). The Pearson's product moment correlation coefficient r is shown in each plot.

Melt pond prediction model

As a first step in the development of a f _p prediction model, we considered the importance of evaluating data collected in different conditions separately in this study. A robust f _p prediction model, i.e. one that includes a wide range of FYI and MYI features, would include input data of all features captured in the same controlled conditions. Unfortunately the limiting factor in this study is the scarcity of high-resolution optical scenes from which to derive spring f _p data coincident to winter C-band SAR imagery. Moreover, C-band σ _hh^° and σ _hv^° are affected by θ; both parameters generally decrease with increasing θ, with the rate of decrease greater for FYI compared with MYI (Geldsetzer and others, Reference Geldsetzer, Mead, Yackel, Scharien and Howell2007). This effect would mask backscatter to f _p relationship if not controlled.

Based on the above, we fit linear regression models to the Sentinel-1 backscatter and texture data, and f _p data from the VN site only. Our rationale for choosing data from VN only is based on the ideal timing of those VIS-NIR acquisitions relative to the sub-stage of melt pond evolution. Based on in situ studies (Eicken and others, Reference Eicken, Krouse, Kadko and Perovich2002; Polashenski and others, Reference Polashenski, Perovich and Courville2012 and Landy and others, Reference Landy, Ehn, Shields and Barber2014), rapid changes in f _p linked to melt rates and significant meltwater inputs from snow are known to occur. After the snow cover has ablated, f _p is comparatively stable in time and varies in space according to topographic relief (macroscopic flaws such as seal holes and leads excepted). Eventually permeability thresholds are crossed and rapid vertical drainage of melt ponds occur in a spatially heterogeneous manner. Based on the timing of acquisition after the onset of ponds (14 d), visual inspection of the pan-sharpened VIS-NIR imagery, and statistics derived from classified VIS-NIR scenes (Fig. 8), we are confident that site VN is most representative of this sub-stage of melt pond evolution, when f _p is linked to ice topography, hence closely related to ice type. On the other hand, VIS-NIR data at AB and VS were acquired during the initial snow melt and drainage sub-stages, respectively. The VN data capture winter FYI and MYI that evolves into a wide range of f _p, from 0.07 to 0.73 (see Fig. 8). Regression model fitting was done using two approaches: (i) using σ _hh^° only; and (ii) using σ _hh^° and candidate texture parameters from Figure 10. A table of the output regression model parameters is given in Table 6 and the derived regression equations are shown below.

(1)$$f_{\rm p} = - 0.317 - 0.039X_1,$$

(2)$$f_{\rm p} = - 0.533 + 0.853X_1 - 1.157X_2 - 0.069X_3.$$

Table 6. Summary of regression model outputs. Shown for each model are the coefficient of determination (r ²), standard error of the estimate (S), the F-value, the P-value, and the input predictor variables

The Sentinel-1 data used for regression model development was collected over a θ-range of 36.5–39.7°, a portion of the studied angles within which strong σ _hh^° and texture relationships with f _p were observed. We expect the models would be applicable to a wider range of θ, as the change in σ _hh^° from FYI and MYI is not significant from ~35 to 45° (testing is required). Furthermore, the models may prove capable of predicting f _p of DFYI to within an acceptable error, as both MYI and DFYI display similar σ _hh^° and texture characteristics, and evolve to overlapping f _p in spring. Fortunately, with Sentinel-1 operating as a two satellite constellation with near-polar orbits 180° apart, the θ limitation of the models presented here can be overcome by the regular collection of scenes during the winter period.