1. INTRODUCTION
Glaciers and ice caps (GIC) around the world are changing rapidly, especially in the current global warming context (Vaughan and others, Reference Vaughan, Stocker, Qin, Plattner, Tignor, Allen, Boschung, Nauels, Xia, Bex and Midgley2013). Monitoring them is of primary importance to comprehend and predict the impacts of these changes (e.g. sea-level rise and water resources). Large-scale monitoring of GIC elevation changes and mass balances using geodetic methods has improved in recent decades thanks to new remote sensing techniques that allow the generation of DEMs. At the forefront, optical satellite missions (i.e. ASTER, SPOT-5, ICESat and more recently Pléiades and WorldView) have been extensively used for glaciological purposes (Kääb, Reference Kääb2008; Berthier and others, Reference Berthier, Schiefer, Clarke, Menounos and Rémy2010; Gardner and others, Reference Gardner, Moholdt, Arendt and Wouters2012; Bolch and others, Reference Bolch2013; Trüssel and others, Reference Trüssel, Motyka, Truffer and Larsen2013; Quincey and others, Reference Quincey2014; Papasodoro and others, Reference Papasodoro, Berthier, Royer, Zdanowicz and Langlois2015). However, optical technologies are dependent on weather and lighting conditions, which complicate their systematic and/or on-demand use (i.e. at precise dates or years) in regularly cloudy areas. Given the all-weather capabilities of the C and X bands (i.e. penetration through rain and clouds), there is a need to explore in detail the data derived from Synthetic Aperture Radar (SAR) satellite missions for glaciological measurements.
The two main approaches to generate a DEM from SAR data are interferometry and stereo radargrammetry (SRG). Interferometry uses the phase signal and has long been prioritized given that most spaceborne radar sensors over the past few decades (e.g. ERS-1 and RADARSAT-1) had low spatial resolutions compared with their high resolution in phase (Rott, Reference Rott2009; Capaldo and others, Reference Capaldo, Crespi, Fratarcangeli, Nascetti and Pieralice2011). SRG uses the amplitude signal and has been increasingly used following the launch in recent years of very high-resolution SAR satellites, such as TerraSAR-X, COSMO-SkyMed and RADARSAT-2 (Toutin and others, Reference Toutin, Chenier, Schmitt and Zakharov2009; Raggam and others, Reference Raggam, Gutjahr, Perko and Schardt2010; Toutin, Reference Toutin2010, Reference Toutin2012; Capaldo and others, Reference Capaldo, Crespi, Fratarcangeli, Nascetti and Pieralice2011, Reference Capaldo2015; Toutin and Omari, Reference Toutin and Omari2011). These recent SRG studies have shown promising results with vertical precisions of up to ~4 m (standard deviation, SD) for elevations extracted over flat areas. This precision is somewhat worse than the relative error of the interferometry-derived global TanDEM-X DEMs (i.e. better than 2 m; see Krieger and others, Reference Krieger2007; Bräutigam and others, Reference Bräutigam2015). However, the slightly weaker precision of SRG is somewhat compensated for by its higher operability and capability to quickly generate on-demand DEMs with no ground control point (GCP), thanks to the Rational Polynomial Coefficients (RPCs), orbit information and accurate geometric models. By contrast, along with the fact that TanDEM-X DEMs will not be available for on-demand dates, the interferometric approach for DEM generation can be affected by a lack of coherence and decorrelation if the SAR images are not acquired simultaneously.
An advantage of using RADARSAT-2 (hereafter noted R2), compared with COSMO-SkyMed or TerraSAR-X, is that each R2 image is systematically provided with RPCs. This allows direct and fast geometric modeling without GCP in PCI Geomatics (2013)© software. For both TerraSAR-X and COSMO-SkyMed, RPCs are not distributed with images and a tool for RPC generation has to be used after acquisition (using for example the Italian SISAR software; Capaldo and others, Reference Capaldo, Crespi, Fratarcangeli, Nascetti and Pieralice2011, Reference Capaldo2015).
Given the above advantages, R2 SRG could be a good alternative in remote and frequently cloudy regions. This was shown during the finalization of the Canadian Arctic 1:50k map, north of 81°N, in which R2 SRG has been a very efficient replacement for typical cartography methods, such as aerial photogrammetry and optical stereoscopic satellites (Clavet and others, Reference Clavet, Toutin and Kharbouche2011). In this project, the temporal resolution of image acquisition was reduced due to the convergence of the polar sun-synchronous orbits of R2 at very high latitudes (Clavet and others, Reference Clavet, Toutin and Kharbouche2011; Toutin and others, Reference Toutin, Blondel, Clavet and Schmitt2013). On the other hand, the R2 SRG in C-band (5.4 GHz) can be affected by different limitations, notably a varying penetration depth on ice and snow depending on the date and surface conditions (Dall and others, Reference Dall, Madsen, Keller and Forsberg2001; Rignot and others, Reference Rignot, Echelmeyer and Krabill2001; Kääb and others, Reference Kääb2015), geometric distortions due to the side-looking view of the SAR acquisitions (shadow, layover and foreshortening), and radiometric complications (i.e. lack of texture and speckle noise). In very challenging mountainous areas, these limitations could decrease the accuracy, vertical precision and efficiency of R2 SRG derived products (Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012, Reference Toutin, Blondel, Clavet and Schmitt2013). Nevertheless, as the present paper shows, the impact of these constraints can be limited when acquisitions are made at an optimal time of year, with ideal acquisition geometry and for low-to-medium relief.
Here, the potential of the R2 SRG (with the Wide Ultra-Fine (WUF) image mode) is presented for DEM extraction, as well as for elevation changes calculation (historical and recent) on Barnes Ice Cap (BIC) (Nunavut, Canada). To our knowledge, this is the first use of R2 SRG exclusively for glaciological purposes. A comparison between two geometric models, namely hybrid Toutin's model (HTM) and rational function model (RFM), is also carried out to determine which is the most accurate for geometric modeling of R2 WUF images without GCP. Finally, an analysis of possible biases (e.g. C-band penetration depth) in R2 SRG-derived DEMs is carried out.
2. STUDY REGION AND DATASETS
2.1 Study region
BIC (70.08°N; 73.66°W) is located in the middle of Baffin Island, the southernmost island of the Canadian Arctic Archipelago (Fig. 1). In 2002, it covered an area of ~5860 km2 (Pfeffer and others, Reference Pfeffer2014), with a summit reaching 1100 m a.s.l. and the lowest elevations at <400 m a.s.l. BIC is a flat ice cap with an average slope of ~1.6°, while the mean slope of the surrounding ice-free terrain is ~3.6°. It is worth mentioning that BIC accumulates exclusively by superimposed ice (Zdanowicz and others, Reference Zdanowicz2012; Dupont and others, Reference Dupont2012; Reference Dupont2014). Previous studies and winter field works conducted near the summit in March 2011 by Dupont and others (Reference Dupont2014) revealed no firn and a seasonal snowpack ~1 m thick only. BIC, along with Penny Ice Cap, are partial remnants of the Laurentide Ice Sheet (Zdanowicz and others, Reference Zdanowicz2012). Given the remoteness and the simple geometry of this ice cap, as well as the strong cloudiness of this Arctic region (Orvig, Reference Orvig1954; Svoboda and Paul, Reference Svoboda and Paul2009), BIC represents a relevant case study for R2 SRG. Furthermore, BIC experienced an accelerated thinning in recent decades (see Sneed and others, Reference Sneed, Hooke and Hamilton2008; Dupont and others, Reference Dupont2012; Gardner and others, Reference Gardner, Moholdt, Arendt and Wouters2012) suggesting the possibility of measuring elevation changes over a small time interval.
Fig. 1. Study area superposed by different datasets of the project. The background raster image is a hillshade derived from the CDED.
2.2 R2 images
The ten R2 images used in this project (i.e. five stereoscopic pairs, Fig. 1; Table 1) were acquired through the Canadian Space Agency's SOAR program (Science and Operational Applications Research). This program allows Canadian researchers to freely plan and acquire R2 images for various applications.
Table 1. R2 images and stereoscopic pair characteristics
The acquisition geometry is of prime importance when using SAR images for SRG purposes. To maximize the sensitivity to topography and thus the accuracy of the extracted parallax, the intersection angle between the two acquisitions should be large (Toutin and Gray, Reference Toutin and Gray2000; Méric and others, Reference Méric, Fayard and Pottier2009). However, in the radar spectrum, a large intersection angle between two acquisitions decreases the similarity of the two images and thus complicates the image matching process. The accepted compromise is to use images acquired on the same side with a base-to-height ratio (B/H), which is an indicator of the sensitivity to topography, between 0.35 and 0.70 (Capaldo and others, Reference Capaldo2015). In previous studies that used R2 SRG in regions where there are GIC (Clavet and others, Reference Clavet, Toutin and Kharbouche2011; Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012), shallow angles of ~37–39° and ~48–49° were preferred, which corresponds to B/H between 0.31 and 0.37. This range of incidence angles represents a coherent geometry for minimizing the impact of ice and snow on the C-band signal (e.g. rapid changes affecting radiometry and penetration depth) and thus, the resulting image and DEM processing (Toutin and others, Reference Toutin, Blondel, Clavet and Schmitt2013). Hence, we used similar acquisition geometries for 4 pairs out of 5, namely shallow angles of ~36–37° and ~47°, corresponding to B/H ranging between 0.33 and 0.37 (Table 1, DEMs 1–4). Because of R2 user conflicts, the images encompassing the southern dome of BIC were acquired with slightly steeper angles (32.1 and 43.5°, B/H of 0.32). Images that form a stereoscopic pair also had to be acquired on two closely spaced dates (i.e. maximum of 11 d here) to have the surface conditions as stable as possible between acquisitions.
Images were acquired at the end of the 2013 ablation season (i.e. September/October and beginning of November for DEM 5; Table 1 for the exact dates). Thus, the surface texture over the ice cap surface for the images of DEMs 1–4 is pronounced (Fig. 2) and must lead to good correlation matching and DEM coverage. The high-water content of the surface layer at this time of the year is also expected to minimize the C-band penetration depth (Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012, Reference Toutin, Blondel, Clavet and Schmitt2013). These assumptions are somewhat more uncertain for the images of DEM 5, especially over the accumulation area where a fresh layer of snow could probably be present. No other dataset was retrievable to confirm this but the very weak coverage of DEM 5 (Section 4.1) over this area suggests this situation. However, based on the winter surface conditions observed by Dupont and others (Reference Dupont2014), this fresh snow layer is probably very thin at the beginning of November. Moreover, the overlapping between DEMs 4 and 5 (Fig. 1) and the poor coverage of DEM 5 over the accumulation area leads to the fact that, for the penetration depth analysis presented in this paper, this area consists almost exclusively of elevations from DEM 4. Thus, the impact of the surface state of the images of DEM 5 on our penetration depth results (Section 4.2) is negligible.
Fig. 2. (a) Portion of the R2 image acquired on 28 September 2013 (DEM 1); (b) Portion of the R2 image acquired on 23 October 2013 (DEM 3); (c) Portion of a R2 image that was not used in the project (i.e. missing another image to form a stereoscopic pair due to user conflict) and acquired on 17 August 2013. Ablation areas are identified as ‘i’, ice-free terrain is identified as ‘ii’ and accumulation areas as ‘iii’.
All images were acquired in WUF mode (pixel spacing of 1.6 m × 2.8 m (range × azimuth); swath coverage of 50 km × 50 km). Individual scenes acquired using this mode cover an area 6 times larger than Ultra-Fine image mode scenes (20 km × 20 km), which were used in previous assessments of R2 SRG (e.g. Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012, Reference Toutin, Blondel, Clavet and Schmitt2013; Capaldo and others, Reference Capaldo2015). The WUF mode was thus preferred to reduce the number of acquisitions necessary to cover the large BIC. Following Toutin (Reference Toutin2010), the single look complex (slant-range georeference) format and HH polarization were preferred. Both descending and ascending orbits were chosen.
2.3 Historical canadian digital elevation data (CDED)
Historical CDED (scale of 1:50k) covering the ice cap and the surrounding ice-free terrain were retrieved from the Natural Resources Canada web portal (http://www.geogratis.ca/). They were used for historical mass-balance calculations. Derived by stereo-compilation of aerial photos acquired during late summer between 1958 and 1961, CDED elevations are orthometric and referenced to the Canadian Gravitational Vertical Model of 1928 (CGVD1928). The average CDED elevation differences and their SD off glacier between 340 CDED tiles covering Baffin Island and ICESat laser altimetry points were reported to be 1.1 and 5.1 m, respectively (Gardner and others, Reference Gardner, Moholdt, Arendt and Wouters2012). A careful visual interpretation was performed on the shaded relief images derived from the CDEDs covering BIC (not shown here) and suggested that these DEMs have a very good reliability. Meltwater channels are well captured in the topography and no obvious artefact due to interpolation in the accumulation area was observed. Before merging every CDED, each individual tile was projected in UTM projection (WGS84 datum) and resampled at a spatial resolution of 30 m.
2.4 Airborne laser altimetry: IceBridge airborne topographic mapper (ATM)
Airborne laser altimetry data acquired from the ATM NASA IceBridge mission were retrieved for May/June 2005 and 2011. The 2005 data were used for recent mass-balance calculations, while the 2011 data were used for absolute co-registration of the SRG DEMs on ice-free terrain. The laser from the ATM system scans the surface in an approximately conical pattern (Krabill and others, Reference Krabill2002). For the 2005 acquisition system, each laser shot has a footprint of 1–3 m and a spacing of 2–5 m between each shot, whereas the 2011 system uses a laser footprint of ~0.5 m with a measurement density of 1 point per 10 m2 (Gardner and others, Reference Gardner, Moholdt, Arendt and Wouters2012). The 2005 acquisition system allows a nominal accuracy of <0.2 m, while for 2011 it is <0.1 m. Given the flat slope of BIC, the resampled Icessn product was chosen (Gardner and others, Reference Gardner, Moholdt, Arendt and Wouters2012). Data were projected to UTM map projection (WGS84 datum) and converted to orthometric reference according to the Canadian Gravitational Model of 2005 (CGG2005).
2.5 SPOT-5 high-resolution stereo (HRS) DEM
A SPOT-5 HRS DEM derived from images acquired on 31 August 2010 was acquired from the SPIRIT IPY project (Korona and others, Reference Korona, Berthier, Bernard, Rémy and Thouvenot2009). Version 1, which was generated using correlation parameters adapted for gentle terrain, was preferred. Evaluated by Gardner and others (Reference Gardner, Moholdt, Arendt and Wouters2012), the vertical precision of this DEM on stable terrain is 4.9 m (SD). Elevations are orthometric and referenced to the EGM96 geoid. The DEM was projected to the UTM map projection (WGS84 datum) and resampled to a spatial resolution of 30 m. This dataset was used for the penetration depth analysis (Sections 3.2.3 and 4.2).
2.6 Delineation of ice and ice-free terrain
A precise distinction between ice-free and ice-covered terrain is of prime importance when assessing DEMs. To encompass the large BIC and the surrounding terrain, two cloud-free Landsat-5 images (spatial resolution of 30 m) were retrieved from the USGS web portal (http://glovis.usgs.gov/). The image covering the whole BIC and the surrounding southern area is from 10 August 2010, while the image covering only the surrounding northern area is from 5 July 2003. Both images were classified (supervised maximum likelihood) to extract the recent ice-free and ice-covered areas for the whole study zone.
Furthermore, computing glacier mass changes using the geodetic method requires delimitation of the glacier for years that are close to those of the DEMs used (i.e. 1960, 2005, 2013). For the older delimitation of BIC, we used the raw vector files from the 1:250k Canadian National Topographic Data Base acquired in 1958/1960. The 2005 margin was manually digitized from a Landsat-5 image from 27 August 2005 retrieved from the USGS web portal. The most recent BIC limit was taken from the classification made on the Landsat-5 image acquired on 10 August 2010.
3. METHODS
3.1 R2 DEM generation without GCP
All R2 image and DEM processing was performed with the PCI Geomatics (2013)© commercial software. Every R2 image was first filtered in the Focus module using a Frost filter with a 9 × 9 window (Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012), to eliminate speckle. This 9 × 9 filter was chosen following previous tests (not shown here) which were intended to select the filter offering the best compromise between speckle elimination and texture conservation.
SAR images are affected by various geometric distortions related to the platform, the sensor and the Earth for which correction and modeling with geometric models are required (Toutin, Reference Toutin2011). Two different geometric models requiring no GCP were tested in the OrthoEngine module: the HTM and the RFM. The HTM is a recently developed hybrid model. It combines the well-known Toutin's rigorous physical model (see Toutin and others, Reference Toutin, Chenier, Schmitt and Zakharov2009), which uses metadata to approximate parameters of the model, with the empirical RFM (i.e. related to the supplied RPCs), which generates virtual points that are then used as inputs to the model. Further details on the HTM can be found in Toutin and Omari (Reference Toutin and Omari2011). The RFM is a purely empirical model based on RPCs. Essentially, the RFM relates image pixel coordinates (column and line) to terrain coordinates in the form of polynomial function ratios solved with RPCs (PCI Geomatics, 2013). More details on the RFM can be found in Zhang and others (Reference Zhang, Fei, Li, Zhu and Li2010, Reference Zhang, He, Balz, Wei and Liao2011). Both models were used for each stereoscopic pair in order to determine which was most suitable for geometric modeling of the R2 WUF image mode without GCP.
The following steps were performed for each DEM generation. R2 images were converted into quasi-epipolar images with a downsampling factor of 2 (Toutin and others, Reference Toutin, Blondel, Clavet and Schmitt2013; Capaldo and others, Reference Capaldo2015). This reduces the effect of the remaining speckle on the ice-cap surface (Toutin and others, Reference Toutin, Blondel, Clavet and Schmitt2013). The correlation matching was then performed on the quasi-epipolar images using the default OrthoEngine procedure, namely a hierarchical grey-level image matching (mean normalized cross-correlation method with sub-pixel computation of the maximum of the correlation coefficient). As suggested by Ostrowski and Cheng (Reference Ostrowski and Cheng2000), a low level of detail was selected for the DEM generation to increase the correlation matching. The terrain type setting, which did not have a significant impact on our DEMs, was simply set to ‘hilly’. Furthermore, no interpolation was performed to fill DEM data voids, since it may produce false elevations with larger errors (Berthier and others, Reference Berthier2014). DEMs were finally geocoded with 30 m grid spacing.
As regards the altimetric reference system, note that when the HTM is used in OrthoEngine, DEMs are automatically referenced to EGM2008. By contrast, when using the RFM, DEMs are referenced to the WGS84 ellipsoid and have to be converted to EGM2008. Over our study area, average differences between the various geoids (EGM2008, CGG2005, EGM96, CGVD1928) are below 10 cm and were thus neglected.
3.2 R2 DEM evaluation and correction
3.2.1 DEM assessment
DEM accuracy and vertical precision (i.e. biases and SDs) were evaluated on ice-free terrain by comparing with up to 15 000 ATM elevation points. An evaluation was also performed on the ice cap using a similar amount of ATM points to better compare the efficiency and consistency of the two geometric models. Consequently, the evaluation over the ice cap does not represent the real DEM accuracy (and vertical precision) since it includes errors related to both the variability of the possible penetration depth and elevation changes over the period 2011–13. For the sake of comparability with Toutin and others (Reference Toutin, Omari, Blondel, Clavet and Schmitt2012, Reference Toutin, Blondel, Clavet and Schmitt2013), this assessment was done on the raw DEMs (i.e. before co-registration) and discarding absolute elevation differences >100 m. Note that the assessment on the ice cap was carried out independently for the ablation and accumulation areas that were previously manually digitized using the R2 images. The two areas were clearly distinguishable as the ablation area is darker due to low backscattering, while the accumulation area appears much brighter.
3.2.2 Co-registration and bias analysis over ice-free terrain
Given the higher accuracy of the RFM DEMs (see Section 4), only these DEMs were used for the glaciological measurements and are thus analyzed in this section. First, these DEMs were co-registered (X, Y, Z) to the ATM points from 2011 over ice-free terrain using the relationship between the aspect, slope and elevation differences (Nuth and Kääb, Reference Nuth and Kääb2011). Horizontal shifts were not systematic for each DEM. For DEMs 1–4, most of the horizontal offsets ranged between ~1 and 2 DEM pixels (23–55 m), while DEM 5 was better located with easting and northing shifts of −12.2 and 3.8 m, respectively (Table 2). The vertical biases ranged from −2.3 to −7.3 m. The same co-registration approach was applied to the merged CDED for the historical mass-balance estimate.
Table 2. Horizontal and vertical shifts of the SRG DEMs (RFM model) measured and corrected during the absolute 3-D co-registration
Remote sensing derived DEMs can be affected by various biases (e.g. elevation-dependent bias, cross/along track bias) that can seriously alter the glaciological measurements and their interpretation (Berthier and others, Reference Berthier, Arnaud, Vincent and Rémy2006; Nuth and Kääb, Reference Nuth and Kääb2011). For our 5 SRG DEMs, the ice-free terrain does not cover the whole elevation range so a rigorous elevation-dependent bias analysis could not be correctly conducted. Future work on this particular bias should be done in more mountainous areas. Cross/along track biases, which are often observed in satellite stereoscopic DEMs due to the acquisition geometry, are generally quantified on stable terrain (Nuth and Kääb, Reference Nuth and Kääb2011). For four of the five SRG DEMs (DEMs 1–4), the proportion of ice-free terrain is too small to appropriately analyze and correct such biases. DEM 5 is the only one for which the stable terrain portion is large enough to evaluate these biases but, for the sake of homogeneity in the DEM processing, no correction was applied. Nonetheless, a sensitivity test (not shown) revealed that a cross/along track bias correction on DEM 5 would have decreased the elevations in this DEM by <1 m over the ice cap. Once again, the possible biases could be evaluated in the future for other study areas with more ice-free terrain distributed over the entire scene.
3.2.3 Radar penetration analysis over the ice cap
Elevation extraction on ice and snow with the C-band (5.4 GHz) can lead to a penetration depth of up to several meters (Dall and others, Reference Dall, Madsen, Keller and Forsberg2001; Rignot and others, Reference Rignot, Echelmeyer and Krabill2001; Müller, Reference Müller2011). Here a first-order estimate of the penetration depth was obtained by analyzing the elevation differences on the ice cap between the 2013 SRG DEMs and the 2010 SPOT-5 HRS DEM (Korona and others, Reference Korona, Berthier, Bernard, Rémy and Thouvenot2009). After absolute co-registration on ice-free terrain between both datasets, we hypothesized that the remaining elevation differences on the ice cap mostly results from (1) elevation changes between 2010 and 2013 and (2) the C-band penetration depth. To deduce the penetration depth, we corrected the height change (−2.55 ± 0.5 m) measured by Gray and others (Reference Gray2015) between August 2010 (acquisition date of the SPOT-5 HRS images) and October 2013 (average acquisition date of the R2 images). In detail, we adjusted the SPOT-5 HRS DEM by distributing this mean height change as a function of elevation on the ice cap, using the elevation change gradient with altitude between ATM 2005 and ATM 2011. We thus assumed no temporal differences between the five SRG DEMs, given the similar look and texture of all the R2 images used, and we performed this analysis with the merged SRG DEMs. This was done independently for the accumulation and the ablation areas (Fig. 3 for the magnitude of the corrections over the two areas).
Fig. 3. Penetration depth analysis performed independently for the ablation area (green) and the accumulation area (in blue) of BIC. The uncorrected biases (solid curves) represent the raw elevation differences for each elevation band between the SRG DEMs (2013) and the SPOT-5 HRS DEM (2010), after absolute co-registration on ice-free terrain. The corrected biases (dashed curves) represent the elevation differences that were adjusted for the height lowering that occurred between 2010 and 2013 using results from Gray and others (Reference Gray2015). Section 3.2.3 shows more detail on the analysis method.
3.3 Elevation changes and mass-balance calculation
Historical elevation changes (dH) of BIC were then obtained by subtracting the co-registered CDED from the SRG DEMs (SRG DEM 2013 – CDED 1960). To fully characterize SRG's potential limits for glaciological measurements, we subtracted the 2005 ATM points from the SRG DEMs (SRG DEM 2013 – ATM 2005) to obtain a recent measurement. Since the SRG DEMs were derived from images acquired at the end of the ablation season and the ATM points are from May, a seasonal correction was applied to the ATM points using a value measured by Gray and others (Reference Gray2015) based on Cryosat-2 altimetry measurements. To proceed, we found the recent summer that was the most similar to the summer of 2005. Temperature anomalies in the summer of 2013 were found to be similar to 2005 (based on the Homogenized Canadian Climate Data at the Dewar Lakes station; Vincent and others, Reference Vincent, Zhang, Bonsal and Hogg2002). Thus, the ice cap wide elevation change measured by Gray and others (Reference Gray2015) between May 2013 and October 2013 was added to the 2005 ATM points after distributing according to the 2005–2011 elevation change gradient.
The steps followed to convert the elevation changes to a mass balance are described elsewhere (Papasodoro and others, Reference Papasodoro, Berthier, Royer, Zdanowicz and Langlois2015) and thus not repeated here. For the volume to mass conversion, a density of 900 kg m−3 was used, to facilitate direct comparison with results from Gardner and others (Reference Gardner, Moholdt, Arendt and Wouters2012). As regards mass-balance accuracies, the CDED values were already shown to be highly autocorrelated (Gardner and others, Reference Gardner, Moholdt, Arendt and Wouters2012), so the dH uncertainty of the historical mass-balance calculation was assumed equal to the SD measured on ice-free terrain between the CDED and the SRG DEMs. As a matter of uniformity, the same conservative approach of dH uncertainty was adopted for the 2005–2013 period. For comparison with results from Gardner and others (Reference Gardner, Moholdt, Arendt and Wouters2012), the same uncertainties were assigned for ice density (±17 kg m−3), ice-cap area (10%) and extrapolation to the missing pixel elevations (10%).
4. RESULTS AND ANALYSIS
4.1 Evaluation of the R2 SRG DEMs
SRG DEMs and validation points are shown in Fig. 1 (only DEMs derived from RFM are shown for the map readability). Over ice-free terrain, more than 85% of the area is extracted, regardless of the geometric model, due to the good texture and high resolution of the R2 images. Over the ice cap, this proportion decreases to ~33% and ~26% for DEMs derived from the RFM and HTM, respectively. In the ablation area, the very pronounced texture and radiometric contrast (Fig. 2), related to a wet ice surface, improves the success of correlation matching and thus increases the covered area (>50%). By contrast, the soaked snow and the superimposed ice over the accumulation area (higher than ~800–850 m a.s.l.) negatively impacts the matching and elevation extraction for DEMs 1–4 (~11% of covered area). For DEM 5, the coverage over the accumulation area is even lower (~4%), a poor coverage likely due to the presence of fresh snow on the ice-cap surface as discussed in Section 2.1.
Average accuracy and vertical precision (i.e. bias and SD, respectively) of the raw SRG DEMs measured over ice-free terrain and over the ice cap are given in Table 3 for both the HTM and the RFM. The vertical precision of the extracted DEMs is 7.3 m on average for the RFM DEMs (i.e. between 5.8 and 7.1 m for DEMs 2–5 and 11 m for DEM 1). The vertical precision is 8.5 m on average for the HTM DEMs (i.e. between 7.4 and 11.7 m for the 5 DEMs). These SD are significantly lower than those obtained over a Canadian Arctic mountainous region (SD of up to 30 m; Toutin and others, Reference Toutin, Blondel, Clavet and Schmitt2013) and just slightly higher than the precisions of 4–5 m (SD) measured on flat areas north of Québec City, Canada (Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012; Capaldo and others, Reference Capaldo2015). Even though this evaluation was done prior to co-registration, it is worth noting that the co-registration of the DEM had no significant impact on the SD (<0.2%). Over the ice cap, the SD measured between the validation points (ATM 2011) and the raw SRG DEMs cannot readily be interpreted in terms of vertical precision since it includes errors due to penetration depth variability and elevation changes between 2011 and 2013. Nonetheless, this additional assessment helps to compare the consistency of the two geometric models. Over the ablation area, the SD increases to ~12–13 m. Taking into account the errors mentioned previously (i.e. possible penetration depth and 2011–2013 height change), these values are consistent with our precisions measured on stable terrain and still represent an improvement to the SD of 17.9 m obtained on a flat glacier area by Toutin and others (Reference Toutin, Blondel, Clavet and Schmitt2013). This improved precision may be due to the interpolation conducted in the study by Toutin and others (Reference Toutin, Blondel, Clavet and Schmitt2013) and to the very pronounced texture of BIC over the ablation area in our R2 images. Furthermore, the average SD increases up to ~24 m over the accumulation area of BIC, likely because of the low texture of the images in this zone, which leads to a poorer correlation. This issue could have been eliminated by having stereoscopic pairs acquired during the ablation season (i.e. August) when the texture was enhanced (Fig. 2c).
Table 3. Validation statistics of the SRG DEMs generated using the rational function model and Toutin's hybrid model over ice-free terrain and over BIC. The validation is made by comparing the 2011 ATM points
The choice of the geometric model has a pronounced impact on the vertical bias of the DEM derived from WUF R2 images. On stable terrain, the average bias measured with the validation points is −4.0 m for the raw RFM DEMs and 5.8 m for the raw HTM DEMs. Moreover, unlike DEMs derived from the RFM, the vertical bias strongly varies for single HTM DEMs (i.e. between 29.1 and −8.3 m for the 5 DEMs) and even within each DEM. A reason that could explain this inaccuracy is the bigger image size used here (i.e. WUF mode, 50 km × 50 km) compared with that for which the HTM was firstly calibrated (i.e. Ultra-Fine mode, 20 km × 20 km; Toutin and others, Reference Toutin, Chenier, Schmitt and Zakharov2009; Toutin and Omari, Reference Toutin and Omari2011). The loss of accuracy occurring when using the HTM with WUF R2 images limits the use of this geometric model for glaciological measurements in remote areas where GCPs are difficult to acquire. Thus, these results suggest a better consistency and accuracy of the RFM when using this image mode. The relatively small remaining biases are easy to eliminate if a reference dataset (e.g. ICESat) is available over the surrounding ice-free terrain. Hence, only DEMs derived from this geometric model (RFM) are used for the subsequent analysis (i.e. C-band penetration and glaciological measurements).
4.2 First-order estimate of the C-band penetration depth
A first-order estimate of the possible C-band penetration depth over BIC was obtained using a 2010 SPOT-5 HRS DEM corrected for the 2010–2013 height lowering, with the help of results from Gray and others (Reference Gray2015) (Section 3.2.3). As shown in Fig. 3, the SRG DEM remains slightly below the corrected SPOT-5 DEM over both the ablation area (0.23 ± 0.66 (SD) m below SPOT-5) and the accumulation area (0.07 ± 1.54 (SD) m below SPOT-5). Given the SD of the two values, we consider the estimated penetration depth to be insignificant and not distinguishable from 0. Thus, no correction was applied to the SRG DEMs. These low values are likely related to a high dielectric constant (i.e. humid ice surface on lower parts and soaked snow/superimposed ice at higher elevations) that occurs at the end of the ablation season, attenuating the radar signal and leading to a surface response instead of a volume response (Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012, Reference Toutin, Blondel, Clavet and Schmitt2013). Our results are also consistent with other values in the literature. Over an exposed ice surface in Greenland, Rignot and others (Reference Rignot, Echelmeyer and Krabill2001) found a very weak C-band penetration of 1 ± 2 m. In addition, Dall and others (Reference Dall, Madsen, Keller and Forsberg2001) measured no C-band penetration over a soaked zone in East Greenland.
4.3 Elevation changes and mass balances
Figure 4 shows historical (1960–2013) and recent (2005–2013) elevation changes (dH) for BIC. The historical dH measured for the 5 DEMs range between −26 and −42 m (i.e. annual elevation changes (dH/dt) of −0.5 to −0.8 m a−1). The majority of the areas located below 600 m a.s.l. experienced thinning rates of 1 m a−1. Between 600 and 1000 m a.s.l., an altitude range encompassing ~70% of BIC area, dH/dt is ~−0.5 m a−1. For the recent period from 2005 to 2013, the glacier-wide dH/dt decreased significantly to −1.2 m a−1 with a strong gradient with elevation (Fig. 2c). Between an elevation of 400 and 1100 m a.s.l. (97% of BIC area), thinning rates of 3.6 m a−1 to 0.5 m a−1 are measured. These values are in agreement with Gardner and others (Reference Gardner, Moholdt, Arendt and Wouters2012) and Sneed and others (Reference Sneed, Hooke and Hamilton2008), and confirm the clear recent acceleration of the thinning of BIC.
Fig. 4. (a) Rate of elevation change (dH/dt) during 1960–2013 between CDEDs and SRG DEMs. (b) dH/dt during 2005–2013 between 2005 ATM points and SRG DEMs. (c) dH/dt as a function of altitude for the two study intervals.
Glacier-wide mass balances are −0.52 ± 0.19 m w.e. a−1 (−3.1 Gt a−1) during 1960–2013 and −1.06 ± 0.84 m w.e. a−1 (−6.2 Gt a−1) during 2005–2013 (Table 4), which is a doubling of the specific rate of mass loss. Our mass balances are very similar to Gardner and others (Reference Gardner, Moholdt, Arendt and Wouters2012) for analogous historical and recent periods (Table 4), again confirming the recent strong acceleration of mass loss over BIC. Although our error bars are high, our results reveal the potential of using R2 SRG (WUF images) for relatively short-time intervals when dH is significant.
Table 4. Mass balances of the BIC calculated for the periods 1960–2013 and 2005–2013. On the right side, results of Gardner and others (Reference Gardner, Moholdt, Arendt and Wouters2012) for similar year intervals are shown for comparison
5. CONCLUSIONS AND PERSPECTIVES
Here we explored the potential of SRG for glaciological purposes using very high resolution R2 WUF images (pixel spacing of 1.6 m × 2.8 m; swath coverage of 50 km × 50 km). We assessed the quality of SRG DEMs derived from two geometric models that require no GCP, namely the HTM and the RFM. Over the low relief ice-free terrain surrounding BIC, average vertical precisions of 7.3 and 8.5 m (1σ confidence level) were measured from the RFM and HTM DEMs, respectively. In addition, RFM DEMs appeared more accurate than the HTM DEMs with an average vertical bias of −4 m (i.e. between −2 and −7 m), in comparison with biases ranging between 30 and −8 m for the HTM DEMs. The HTM DEMs also showed spatially varying elevation biases even within a single DEM, thus complicating the corrections and the use of the HTM DEMs for mapping glacier elevation changes. Hence, our analysis suggests that the RFM provides a more appropriate geometric model for R2 WUF images. Our vertical precisions are close to the precisions of 4–5 m obtained in a more optimal study site (Toutin and others, Reference Toutin, Omari, Blondel, Clavet and Schmitt2012; Capaldo and others, Reference Capaldo2015). Our mass-balance estimates have shown an acceleration in mass loss in recent years, which is in agreement with other recent studies (Sneed and others, Reference Sneed, Hooke and Hamilton2008; Gardner and others, Reference Gardner, Moholdt, Arendt and Wouters2012). This demonstrates, for the first time, the efficiency of R2 SRG for glaciological measurements over both long- and short-time intervals when thinning or thickening rates are significant.
Our first-order analysis suggests a negligible C-band penetration depth over BIC in our SRG DEMs derived from R2 images acquired at the very end of the ablation season (late September to beginning of November). This result is in agreement with former studies (Dall and others, Reference Dall, Madsen, Keller and Forsberg2001; Rignot and others, Reference Rignot, Echelmeyer and Krabill2001) and is related to a wet surface (humid ice, soaked snow and superimposed ice) that leads to surface radar response. This is an important result showing the added value of the R2 SRG technology for retrieving elevation changes when images are acquired at a suitable period of the year. In order to analyze the potential of SRG for every glacier accumulation regimes, further tests of SRG should be conducted on glaciers with cold firn persisting throughout the ablation season in the upper reaches.
According to our findings, the acquisition date of the images is the main constraint to focus on for an optimal use of R2 SRG for glaciological purposes. It is of primary importance to acquire images during or directly at the end of the ablation season to maximize the radiometry and texture of the image (Fig. 2), as well as to reduce the penetration depth. Initially, our project involved the use of R2 images acquired exclusively at this time of year (i.e. DEMs 1 and 2; as well as Fig. 2c). However, some of our scheduled acquisitions were impeded due to user conflicts and this led to new image acquisition later in the autumn. The major impact of those belated acquisitions was the less pronounced texture of the images over the ice cap, which decreased the coverage of the derived DEMs. Further tests could be conducted by acquiring R2 images directly during the ablation season (e.g. August for the Canadian Arctic) to analyze the impact in terms of DEM coverage. Moreover, we stress that the difference in the dates between two images of a stereoscopic couple must be as small as possible (e.g. no more than 11 d in the case of R2) to minimize changes in the surface state between the images and to avoid any surface displacement due to ice flow between the two dates. This is of major importance to avoid any undesired parallax that could negatively affect the generated DEMs, especially over fast flowing glaciers. It is expected that the upcoming RADARSAT Constellation Mission (RCM) will reduce the aforementioned limitations. Planned to be completely launched in 2018, this constellation of three identical satellites will offer greater revisit capabilities than R2 alone, especially in the northern regions that will be covered up to four times daily (Chalifoux, Reference Chalifoux2015). Therefore, the possibilities of acquisition conflicts between users of the RADARSAT program will necessarily be reduced by the RCM.
In terms of operability, the R2 SRG technology has considerable advantages for the study of glaciers located in remote and cloudy regions. First, the all-weather capability of the C-band allows image acquisition even in the presence of clouds or fog. This adds to the capability of geometric treatments without the need of GCPs, which are often difficult to obtain in such remote regions. The small remaining bias within DEMs derived from the RFM can be eliminated by co-registering to any good reference dataset. Given the above, GIC in Patagonia and Alaska (and other maritime environments) could be appropriate targets for further applications of R2 SRG for glaciological measurements, given their remoteness and the cloudy conditions that characterize those regions. Since the results of the present paper are based on a relatively flat ice cap that represents an optimal case study, further tests of SRG should be conducted for glaciers located in steeper reliefs to fully characterize the glaciological potential of the SRG.
ACKNOWLEDGEMENTS
Charles Papasodoro acknowledges support from the Fond Québécois de Recherche en Nature et Technologies (FQRNT) fellowship program. We thank the Canadian Space Agency for freely providing the RADARSAT-2 images through the SOAR program. The National Research Council of Canada supported this project. Étienne Berthier acknowledges support from the French Space Agency (CNES) through the TOSCA program. We thank Alex S. Gardner for his feedback that led to an improvement of this work and Christian Zdanowicz for helpful discussions on this research. We also thank Catherine Brown, from the Geomatic Department of the University of Sherbrooke, for improving the language in the manuscript.
Open access
