Area , elevation and mass changes of the two southernmost ice caps of the Canadian Arctic Archipelago between 1952 and 2014

Introduction Conclusions References


Introduction
With a glacierized area of ∼ 150 000 km 2 , the Canadian Arctic Archipelago (CAA) is one of the major glacier regions in the world (Pfeffer et al., 2014).In response to recent Arctic warming (Tingley and Huybers, 2013;Vaughan et al., 2013;Comiso and Hall, 2014), glaciers in the CAA have experienced an acceleration in their mass loss.For the southern parts of the CAA, annual thinning of glaciers has doubled between the historical  and recent (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) periods (Gardner et al., 2012).Over the entire CAA, the rate of mass change has tripled between 2004 and 2009, reaching −92 ± 12 Gt a −1 during the period 2007-2009 (Gardner et al., 2011), making this region one of the main contributors to eustatic sea-level rise for this period, after Greenland and Antarctica (Gardner et al., 2013;Vaughan et al., 2013).Continued monitoring of CAA glaciers is thus critical.
Located in the southeastern part of the CAA, Baffin Island is the largest island of the archipelago (Andrews et al., 2002) and has a total ice-covered area of ∼ 37 000 km 2 .In addition to two major ice caps, Barnes (∼ 5900 km 2 ) and Penny (∼ 6400 km 2 ), Baffin Island is also covered by a num-ber of isolated icefields and small ice caps, including Grinnell Ice Cap (GRIC) and Terra Nivea Ice Cap (TNIC) on Meta Incognita Peninsula, at the southernmost tip of the island (Fig. 1).Compared to Barnes and Penny Ice Caps, GRIC and TNIC have received little scientific attention so far (Andrews et al., 2002).Different in situ geophysical measurements were carried out in the 1950s (Blake, 1953;Mercer, 1956), and in the 1980s by teams from the University of Cambridge and the University of Colorado (Dowdeswell, 1982(Dowdeswell, , 1984)).Under the supervision of veteran Norwegian glaciologist Dr. Gunnar Østrem, other measurements were conducted on GRIC in the early 1990s by a scientific team from Bates College (Maine, USA) and Nunavut Arctic College.Later in [2003][2004], glaciologists from the Geological Survey of Canada carried out in situ measurements on GRIC (Global Navigation Satellite System (GNSS) elevation measurements, automatic weather observations, snow depth and surface mass balance measurements) with the objective of establishing a long-term observing site.However, consistent prohibitive weather conditions coupled with difficult access to the ice cap led to the cancellation of the project (Zdanowicz, 2007).A recent study (Way, 2015) analyzed changing rates of glacier recession for GRIC and TNIC since the 1950s using historical aerial photographs, satellite (Landsat) imagery and digital elevation models (DEMs).In the present study, we supplement these results by presenting a comprehensive analysis of historical and recent fluctuations of area, surface elevation and mass for GRIC and TNIC over the period 1952-2014.This is done by combining data from spaceborne instruments (laser altimetry and optical stereo imagery), DEMs, airborne imagery (air photos) and in situ GPS surveys.Our analysis differs from that of Way (2015) in the choice of photos, DEMs, and spaceborne, remotely sensed data used to determine glacier change.In particular, we explored the use of sub-meter resolution stereo pairs obtained from the Pléiades satellites to derive DEMs and to collect accurate, numerous and homogeneously distributed ground control points (GCPs) for the photogrammetric processing of aerial photos.We place our findings in the context of the observed pattern of regional glacier changes across the CAA, and discuss climatic forcing factors of particular relevance for the southernmost Baffin Island region.

Study area
GRIC and TNIC (Fig. 1) are located on Meta Incognita Peninsula, 200 km south of Iqaluit, Nunavut.GRIC (62.56 • N, 66.79 • W) covers an area of 107 km 2 (August 2014; this study) with the highest elevations reaching 800 m a.s.l.(above sea level).On the northeast side, some outlet glaciers extend to Frobisher Bay, which opens into the Labrador Sea.TNIC (62.27 • N, 66.51 • W) is located ∼ 17 km south of the GRIC.It covers an area of approximately 150 km 2 (August 2014; this study) with a similar el- evation range as GRIC.Mercer (1954) suggested three factors supporting the continued presence of plateau ice caps on Meta Incognita Peninsula: (1) cool summers, (2) frequent low-level cloud and (3) heavy snowfall.Data from the permanent weather station in Iqaluit (34 m a.s.l.) indicate that winter temperatures (DJF) in this region averaged −24 • C over the past 60 years, while mean summer temperatures (JJA) averaged 6.5 • C. Total annual precipitation is ∼ 500 mm (snow: ∼ 300 mm; rain: ∼ 200 mm).Field observations in winter 2003-04 found no firn at the summit of GRIC, and the estimated winter snow accumulation there (∼ 2-3 m snow; or ∼ 0.65-0.75m water equivalent) was approximately equal to the amount of melt in summer (Zdanowicz, 2007).Hence the summit of the GRIC is probably close, or slightly below, the present-day equilibrium line altitude (ELA), making it highly susceptible to experience net mass losses (Pelto, 2010).
Observations from various expeditions in the 1950s revealed that the western margin of GRIC was relatively stable, but that coastal outlet glaciers (eastern margin) were shrinking moderately when compared to photographs from 1897 (Mercer, 1954(Mercer, , 1956)).Moraines studied near both ice caps in the early 1980s indicated that the most recent phase of recession dated from the last 100 years, and that both ice caps The Cryosphere, 9, 1535-1550, 2015 probably reached their largest areal extent during the Little Ice Age cold climate interval (Muller, 1980;Dowdeswell, 1982Dowdeswell, , 1984;;Andrews, 2002).Dowdeswell (1982) estimated that the outlet glacier of GRIC that calves into Watts Bay extended much further out a few centuries earlier, but also reported that another outlet glacier to the south of the ice cap was advancing.

Pléiades stereoscopic images
Launched on 17 December 2011 and 2 December 2012 respectively, the Pléiades 1A and 1B satellites have recently shown their high potential for glacier DEM extraction and thus, for mass balance estimations (Wagnon et al., 2013;Berthier et al., 2014;Marti et al., 2015).The two satellites follow the same near-polar sun-synchronous orbit and provide panchromatic and multispectral imagery at a very high ground spatial resolution, 0.7 m for panchromatic and 2.8 m for multispectral images, respectively (Astrium, 2012).Both satellites have independent stereoscopic capabilities.The fact that the panchromatic band images derived from Pléiades satellites are coded in 12 bits represents a clear advantage on a glacier surface (especially over the low contrast accumulation area), given the fact that a large radiometric range provides better contrast and reduces the risk of image saturation (Berthier et al., 2014).
Three stereoscopic pairs were acquired over our study area (Table 1): one for GRIC (3 August 2014) and two for TNIC (14 August 2014 for the eastern part and 26 August 2014 for the western part, with an overlapping area of 84 km 2 ).The stereoscopic pair covering GRIC is cloud-free, while a few clouds (< 10 % of the scene) were present over TNIC during scene acquisitions (Fig. 2).Acquisitions were made at the end of the ablation season to ensure a maximum degree of surface texture (Berthier and Toutin, 2008).Each image was provided with Rational Polynomial Coefficients (RPCs), which allows geometric modeling without GCP.Stereoscopic pairs were used (1) for DEM generation on both ice caps and (2) for GCP extraction for the photogrammetric processing of the historical aerial photos on GRIC (see Sect. 4.2).

Historic Canadian Digital Elevation Data
Historic Canadian Digital Elevation Data (CDED, Natural Resources Canada), provided at a scale of 1 : 50 k, were acquired for the two ice caps.These elevations were derived by stereo-compilation of aerial photos acquired during the summers of 1958 and 1959.Raw elevations are orthometric and referenced to the Canadian Gravitational Vertical Model of 1928 (CGVD1928).The average elevation differ-ences and their standard deviation (SD) between CDED and ICESat laser altimetry were previously calculated off-glacier for 340 CDED maps tiles covering Baffin Island and were reported to be 1.1 and 5.1 m, respectively (Gardner et al., 2012).Here, CDED were used (1) as historical elevations for TNIC and (2) elevations of the surrounding ice-free terrain were used for absolute co-registration for both ice caps (see Sect. 4.3).Artefacts (unrealistic elevations) located in the accumulation area of TNIC were manually identified and deleted using a shaded relief image derived from the DEM.These artefacts were likely due to the poor contrast and low texture of the 1958/59 aerial photos used to generate the CDED.

Historical aerial photos
Historic aerial photos covering GRIC were obtained through the Canadian National Air Photo Library (Natural Resources Canada).We used 24 photos acquired at the end of the ablation season, on 21 and 22 August 1952.A Williamson Eagle IX Cone 524 camera type with a focal length of 152.15 mm was used and the flight altitude was 16 000 ft (∼ 4877 m).The photos are distributed in 3 parallel flight lines with an overlapping coverage of ∼ 30 % between each line and ∼ 60 % between two photos of a same line.These photos, exceptional in their quality of detail and texture, were used for the extraction of historical elevations on GRIC and were preferred to the CDED for this ice cap.In fact, the CDED covering GRIC contained major artefacts (i.e.much larger than for TNIC) in the large snow-covered accumulation zone where the texture was particularly weak on the 1958 photos and thus, was not suitable for historical elevations of this ice cap.

ICESat altimetric points
Surface elevation profiles (GLA14, Release 634) collected by the Geoscience Laser Altimetry System (GLAS) onboard ICESat were acquired (Zwally et al., 2002).Each laser pulsederived footprint corresponds to field-of-view with a diameter of ∼ 65 m and a spacing of 172 m between each footprint (Schutz et al., 2005).ICESat elevations were converted from their original Topex/Poseidon ellipsoid to the WGS84 ellipsoid using tools provided by the National Snow and Ice Data Center.The entire data set (2003 to 2009) was used for absolute co-registration on ice-free terrain, while the data collected during a few selected dates (Table 1) were used to estimate recent elevation changes and to assess the precision of the ASTER August 2007 DEM (see below).

ASTER DEM
Products derived from the ASTER satellite mission have been widely used for glaciological studies (e.g.Kääb, 2008;Nuth and Kääb, 2011;Das et al., 2014).To estimate the recent mass balance for TNIC, we used a DEM (product AST14DMO) generated from an ASTER stereo pair acquired on 3 August 2007.The DEM was automatically derived from bands 3N (nadir-viewing) and 3B (backwardviewing) that have an intersection angle of 27.6 • , which corresponds to a base-to-height ratio of 0.6 (Fujisada et al., 2005).The raw DEM was provided with a grid spacing of 30 m, and elevations are orthometric to the EGM96 geoid.Using 57 ICESat points from two different time periods, namely a few months before (April 2007) and after (November 2007) the ASTER acquisition, we assessed a vertical precision of 2.5 m (SD) on TNIC for this ASTER DEM.Due to cloud coverage, no suitable ASTER DEM was available for GRIC at the end of the ablation season.

In situ GPS measurements
In April 2004, a team from the Geological Survey of Canada measured three surface elevation profiles at 50 m horizontal intervals using a Trimble ® high-precision Real-Time Kinematic GPS system on the southeast, west and northwest sides of GRIC, and at the front of one of its outlet glaciers (Zdanowicz, 2007).Data acquisition was made using a fixed base station on a geodetic benchmark monument, and GPS positions were subsequently processed with the Canadian Center for Remote Sensing's Precise Point Positioning (PPP) System to obtain an accuracy of a few centimeters.For this paper, those transects were used for recent elevation change calculations.Elevations derived from those GPS measurements are referenced to the GRS80 ellipsoid which can be assumed equal to the WGS84 ellipsoid (sub-mm differences).

Glacier outlines
Various data sets have been used to extract the areal extent of the two ice caps at the end of the ablation season (August/September).For GRIC, three data sets from different dates were used.The 1952 outline was derived manually from the orthorectified historical aerial photos.For 1999, we used the ice cap contour from the Randolph Glacier Inventory (RGI 3.2; Pfeffer et al., 2014), which originates from the Canadian CanVec data set for this region, itself derived from a September 1999 Landsat 7 image.We manually digitized the 2014 margin from the orthorectified panchromatic Pléiades image.For TNIC, outlines were derived for four different dates.We used the raw vectors from the 1 : 250 k Canadian National Topographic Data Base as the 1958/59 boundary.Anomalies were found in the delineation of the 1999 margin from the RGI 3.2 (i.e.off-glacier snow patches erroneously included).As an alternative, we manually digi-tized the ice cap margin using a 30 m resolution Landsat 5 image acquired on August 1998.The August 2007 limit was manually traced from an ASTER orthoimage (15 m resolution) provided with the on-demand AST14DMO product, while the 2014 margin was extracted from the orthorectified panchromatic Pléiades images (East and West).To overcome the cloudiness on the Pléiades orthoimages (∼ 20 % of the ice cap outline obscured by clouds), we used a Landsat 8 panchromatic (15 m of resolution) image also acquired in August 2014.The uncertainty assessment of the outlines is presented in Sect.4.4.3.

Meteorological and sea ice records
To quantify changes in the regional climate of the southern Baffin Island region over the period covered in our study, air temperature records were retrieved from the Adjusted and Homogenized Canadian Historical Climate Data of the Iqaluit weather station for the period 1950-2014 (Vincent et al., 2002).This is the permanent weather station in the eastern Canadian Arctic with the most continuous records, extending back to 1946.In addition, time series of sea ice cover area for Hudson Strait and Davis Strait were obtained from the Canadian Ice Service archives over the 1968-2014 period.

Pléiades DEM generation
The Pléiades DEMs were generated using the OrthoEngine module of Geomatica 2013.No GCP were available for the geometric correction so we relied on the RPCs provided with the images.Adding GCP does not improve the vertical precision of the Pléiades DEM, but can reduce the vertical bias (Berthier et al., 2014).The latter bias can be corrected over ice-free terrain when a good reference data set, such as ICE-Sat, is available (Nuth and Kääb, 2011).
The following steps of DEM extraction were repeated for the three Pléiades stereoscopic pairs.First, we collected 20 tie points (TPs) outside and six on the ice cap.Collecting well-distributed TPs was found to improve the relative orientation between the two images providing increased coverage (Berthier et al., 2014).The following processing parameters were used for DEM extraction: the relief type was set to Mountainous and the DEM detail to Low.No interpolation was performed to fill data gaps.Finally, the DEMs were geocoded with a pixel size of 4 m.
Since the ice-free zones on our Pléiades DEM were not large enough to calculate an elevation accuracy with a sufficient number of ICESat points, we report here the vertical precisions obtained in recent glaciological studies.Wagnon et al. (2013) measured a precision of 1 m (SD) on a glacier surface in Himalaya using Pléiades DEM.Berthier et al. (2014)  1 m (SD), highlighting the consistent precision over glacier surfaces.This precision was shown to be mostly correlated with slope.For the small Ossoue Glacier (French Pyrénées), the precision was slightly lower at 1.8 m (Marti et al., 2015).
A similar vertical precision is expected here.

Aerial photos DEM generation
Photogrammetry is widely used in glaciological studies for reconstructing glacier surface prior to the modern satellite era (Fox et Nuttall, 1997;Barrand et al., 2009).In this study, a 1952 DEM of GRIC was created from historical air photos using OrthoEngine.This software uses a mathematical model compensating for both terrain variations and inherent camera distortions (PCI Geomatics, 2013).A typical photogrammetric procedure was then followed to compute the model, solving the least-square bundle adjustment.
Collecting effective GCPs for photogrammetry in mountainous or polar regions remains one of the main difficulties, especially for archive photos (Barrand et al., 2009).To overcome this difficulty, Pléiades-derived products (DEM and orthoimage) were used to collect GCPs.For each aerial stereoscopic model partially covering the surrounding ice-free terrain, 3 to 7 GCPs were collected outside the ice cap on topographic or geomorphologic structures visible on both the Pléiades orthoimage and the aerial photographs.In order to strengthen the mathematical model, every GCP was collected as stereo GCP (i.e. was identified in all possible aerial photographs).A total of 39 stereo GCPs were collected resulting in 106 GCPs.Also, 6 to 10 widely-dispersed TPs were collected for each aerial stereoscopic model.For the models situated in the middle of the photogrammetric block and covering only the ice cap (no ice-free terrain), only TPs were collected in order to connect them to the photogrammetric block.After the final bundle adjustment, the resulting residual averages of all the GCPs were 2.85 m in X, 2.74 m in Y and 2.68 m in Z. TPs residuals were 1.84 m in X and 2.15 m in Y .The generated global DEM was geocoded with a grid resolution of 10 m and no interpolation of data gaps was performed.
Validation of the resulting DEM (before co-registration) against 76 ICESat points on ice-free terrain showed a mean offset of −3.29 m (DEM above ICESat in elevations) and a SD of 4.96 m. Between the Pléiades DEM and the 1952 DEM (co-registered together, see Sect.4.3), the SD of the elevation difference on ice-free terrain was 13.8 m.In total, 66 % of GRIC area was extracted, with data gaps concentrated at the highest elevations in the texture-less accumulation area.

DEM adjustments and co-registrations
DEM co-registration is of primary importance before performing any DEM-based volume change calculations (Nuth and Kääb, 2011).This 3-D co-registration method uses the relationship between aspect, slope and elevation differences over ice-free terrain (Nuth and Kääb, 2011).The Pléiades images only included a small corridor (∼ 20 km 2 ) of icefree terrain near the ice caps (Fig. 2).This corridor coincides with a limited number of cloud-free ICESat points (less than 100 points sparsely distributed around each ice cap) so that a direct co-registration between the Pléiades DEMs and the ICESat points was not optimal.Instead, the CDED tile encompassing the two ice caps was first co-registered with ∼ 1000 ICESat points over ice-free zones.All other DEMs were then 3-D co-registered to the corrected CDED, independently for each ice cap.Thus, the corrected data sets were all referenced to the WGS84 ellipsoid.To evaluate the consistency of the corrections, the offsets over ice-free terrain between each corrected DEM and the corresponding ICESat points were examined.The offset was below 1.5 m in each case, suggesting that the absolute co-registration was well conducted and that the effect of geoid variations (CGVD1928 and EGM96 vs. WGS84) was negligible.However, one must interpret these results with caution given the sparsely distributed and limited number of ICESat points (less than 100) over moderate to hilly terrain.
The two independently co-registered Pléiades DEMs of TNIC (14 and 26 August) were compared in their overlapping zone of 84 km 2 (Fig. 2).The offset measured over the ice-free terrain was −0.1 ± 2.1 (SD) m, while an average elevation difference of 0.64 ± 2.2 (SD) m was measured over the ice cap, probably due to the thinning between 14 and 26 August.These results confirm the robustness of the 3-D co-registration using the corrected CDED.

Elevation changes along ICESat and GPS tracks
For both Grinnell and Terra Nivea Ice Caps, recent elevation changes were measured between 6 ICESat tracks from different laser overpass periods (autumn 2003 to winter 2007) and the 2014 Pléiades DEMs.For GRIC only, elevation changes were also calculated between the April 2004 in situ GPS transects and the 2014 Pléiades DEM.We did not attempt to compute glacier-wide volume or mass changes from those recent elevation changes measurements since (1) they are sparse and only cover a small fraction of the two ice caps and (2) the GPS and some of the ICESat data were acquired at the end of winter, and limited data were available to apply a seasonal correction.Nevertheless, those recent elevation changes along selected tracks were used to complement the differential DEM analysis described below.

DEM-derived elevation changes and mass balances
The geodetic method was applied in order to calculate glacier-wide elevation and mass balances from the DEMs.
The following steps were performed for each calculation.
The Cryosphere, 9, 1535-1550, 2015 First, the co-registered DEMs were subtracted to obtain maps of elevation changes (dH ) and change rates (dH /dt) after dividing by the interval time (dt).The dH values were binned into 50 m elevation bands and averaged after applying a three sigma filter to exclude outliers (Gardner et al., 2012;Gardelle et al., 2013).Pixels with missing data were replaced with the mean dH of the corresponding elevation band.Total volume change for each ice cap (dV ) was then determined by summing volume changes from all elevation bands (n) as follows: where i corresponds to an elevation band of 50 m, H is the mean elevation change measured at elevation band i and A i is the area of the elevation band.In this calculation, the ice cap hypsometry is based on the 1 : 250 k CDED (1958/59), while the ice cap limit conforms to the year of the oldest DEM used in the calculation.Our own sensitivity tests have shown that the choice of the DEM used has a very low impact on the mass balance calculation (< 0.01 m a −1 w.e.), as was also demonstrated in Gardner et al. (2011).
The area-averaged change in elevation over the entire ice cap (glacier-wide), dH /dt avg , was then calculated as follows: where A is the mean of the initial and final ice cap areas, and t is the time interval between the two DEMs.Note that dH /dt (elevation change rate from co-registered DEM subtraction) must be distinguished from dH /dt avg .Finally, the area-averaged specific geodetic mass balance rate (m a −1 w.e.), dM/dt, was calculated as follows: where ρ is the firn and/or ice density.For the historical mass balance of both ice caps, we used ρ = 850 kg m −3 (Huss, 2013), while we used ρ = 900 kg m −3 for the recent period on TNIC (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014).The former value of ρ was chosen assuming that there was a firn zone on the ice cap during the last 6 decades, while a visual interpretation of our images (not shown here) suggests the absence of a significant firn zone after 2007.For the sake of readability, the area-averaged specific geodetic mass balance rate (Cogley et al., 2011) is hereafter simply referred to as mass balance or glacier-wide mass balance.

Accuracy assessment
The main sources of uncertainty in our mass balance estimates are related to uncertainties in the elevation change measurements, the ice cap limits and the density used to convert volume to mass changes.For historical measurements, elevation change uncertainty was assumed equal to the standard deviation over stable terrain between the two coregistered DEMs (GRIC 1952(GRIC -2014: 13.8 m;: 13.8 m;TNIC 1958TNIC -2007: 9.6 m;: 9.6 m;TNIC 1958TNIC -2014: 9 m): 9 m), assuming that elevation errors were 100 % correlated.This is a conservative approach that takes into account both the highly correlated CDED elevation errors (Gardner et al., 2012) and the possible errors associated to the aerial photos-derived DEM (i.e.artefacts and low coverage at higher altitudes).Spatial autocorrelation between the ASTER 2007 and Pléiades 2014 DEMs was analyzed on ice-free terrain to better characterize the elevation change uncertainty in the recent mass balance estimation on TNIC.A low autocorrelation distance (< 100 m) was found between the two elevation products.Applying standard principles of error propagation (e.g.Zemp et al., 2013), we found a very low (±0.1 m) uncertainty for the elevation change averaged over the entire ice cap.Because we consider this value to be likely underestimated, we conservatively assumed that the uncertainty for the elevation changes was equal to the quadratic sum of the two DEMs uncertainties (±1 m for the Pléiades DEM from Berthier et al. (2014) and ±2.5 m for the ASTER DEM from its comparison to ICESat nearly simultaneous on the ice cap), assuming that (1) the elevation errors are fully correlated within each DEM and that (2) errors of the two DEMs are independent.
For ice caps outlines of 1998 and later, we estimated an error of 3 %.This estimate includes possible image interpretation errors (< 2 % of each ice cap extent) and the impact of the image resolution used for outline delimitation (< 1 % of each ice cap extent).Since the ice caps were not covered by debris, interpretation errors were mainly related to the presence of snow-covered surfaces (i.e.snow patches) around each ice cap that may be misinterpreted as ice.The error attributed to the image resolution was established from a comparison analysis made between the Pléiades and Landsat 8-derived TNIC outlines, for which a small difference in extent was obtained (< 1 %).The total uncertainty of 3 % was used for mass balance estimation as well as for area change analysis.For the older (pre-1998) ice cap outlines, a more conservative error of 5 % was applied to take into account the more difficult image interpretation between ice cap limits and snow patches.Those uncertainties are comparable to those reported in Paul et al. (2013) and Winsvold et al. (2014).Finally, an uncertainty of ±60 kg m −3 (Huss, 2013) was assigned to the density factor when estimating the historical mass balance on both ice caps and of ±17 kg m −3 (Gardner et al., 2012) for the recent estimation on TNIC.

Area changes
Areal changes measured for Grinnell and Terra Nivea ice caps since the 1950s are shown in Fig. 2. GRIC experienced a mean rate of areal change of −0.10 ± 0.01 km 2 a −1 between 1952 and 1999, whereas a mean rate of −0.59 ± 0.03 km 2 a −1 is measured for TNIC between 1958/59 and 1998.These rates of area change have been significantly more negative since 1998/99 reaching −1.37 ± 0.04 km 2 a −1 for GRIC and −1.69 ± 0.05 km 2 a −1 for TNIC.For GRIC, the 2014 areal extent is about 20 % smaller than in 1952, while TNIC area shrank by 34 % between 1958/59 and 2014.

Elevation changes
Maps of historical and recent elevation change rates (dH /dt) for the two ice caps are presented in Figs. 3 to 5 (whole ice cap) and Figs. 4 and 6 (changes by elevation).
The glacier-wide rates of elevation change (dH /dt avg ) over the period 1952-2014 were −0.44 ± 0.25 m a −1 for GRIC and −0.56 ± 0.19 m a −1 for TNIC.Similar patterns of historical dH /dt are observed for both ice caps (Fig. 6 and dH /dt maps), revealing an average surface lowering reaching −1.1 ± 0.25 m a −1 for GRIC and −0.9 ± 0.19 m a −1 for TNIC in the lower altitudes (i.e. the outlet glaciers in the peripheries).Above 250 m a.s.l., the thinning rate was consistently 0.1 m a −1 more negative for TNIC than GRIC.The surface thinning was similar for all outlet glaciers of GRIC between 1952 and 2014 (Fig. 3), while on TNIC, a stronger lowering (i.e.< −1 m a −1 ) was experienced on the northeast outlet glaciers between 1958/59-2007 (Fig. 5a).
Elevation change rates sharply increased in recent years for both ice caps.On TNIC, the recent ( 2007 On GRIC, changes in dH /dt were evaluated over the periods 1952-2004 and 2004-2014 using overlapping areas of the 1952 DEM, in situ GPS measurements and ICESat transects ( 2004), and the Pléiades DEM (2014) (Fig. 4).Because of the lack of data about seasonal surface height fluctuations, no correction was applied to account for the different elevation acquisition periods of 1952 and 2014 (August) and 2004 (March/April).For the 203 points where elevation measurements are available for the 3 years (points superposed with black dots on Fig. 4), the dH /dt was up to 6 times more neg- ∼ 11 m for GRIC vs. ∼ 16 m for TNIC) is likely to be at least partly explained by the fact that ICESat transects covering GRIC were located at higher altitudes (mean: ∼ 65 m higher) than those over TNIC.

Mass balances
Mass balances for both ice caps are summarized in Table 2. Between 1952 and 2014, a mass balance of −0.37 ± 0.21 m a −1 w.e. was estimated for GRIC.For TNIC, the historical mass balance (1958/59-2014) was more negative at −0.47 ± 0.16 m a −1 w.e.The mass loss rate for the period 2007-2014 was 5.9 times greater (mass balance: −1.77 ± 0.36 m a −1 w.e.) than that for the period 1958/59-2007 (mass balance: −0.30 ± 0.19 m a −1 w.e.).As previously discussed (Sect.5.2), GRIC has likely experienced a similar acceleration of its mass loss since 2004.6 Discussion

Pléiades as a tool for photogrammetric GCPs collection
In many regions of the world, vast archives of historical aerial photographs represent a potential gold mine for glaciologists in order to document multi-decadal volumetric change of glaciers and ice caps (e.g.Soruco et al., 2009;Zemp et al., 2010).DEMs generated from these aerial photographs allows to put the recent glacier variations (satellite era) in a longerterm perspective.However, these data remain difficult to exploit due to the logistical difficulties involved in the field collection of accurate and well-distributed GCPs in the remote high latitude/altitude regions.Field GCPs were also lacking for the two ice caps studied here.Instead, we took advantage of the very high resolution of the Pléiades imagery (0.7 m) and the vertical precision of the derived DEMs (∼ 1 m) to collect a sufficient number of GCPs for the adjustment of the stereo-model.GCPs were collected on well-defined features that were clearly identifiable on both the Pléiades imagery and the old aerial photos (Fig. 8, yellow arrow).GCP residuals of ∼ 2-3 m in average were obtained after the block bundle adjustment, and a DEM vertical precision of ∼ 5 m (SD) was measured with a few ICESat points available over icefree terrain.This is a satisfactory result given that the aerial photos used here were affected by film distortions that could not be corrected.Our results demonstrate the usefulness of Pléiades stereo-images and DEMs to collect accurate GCPs for photogrammetric processing of old aerial photographs in remote environments.Furthermore, the very fine resolution of Pléiades can help to improve the accuracy of nunataks and/or whole ice caps delimitation, especially when compared to the frequently used Landsat images (Fig. 8).

Comparison to other studies
Our estimates of shrinkage for GRIC and TNIC can be compared with other studies from Baffin Island to verify the coherence of results and get a more complete picture of the pattern of glacier changes across this vast region.Sharp et al. (2014) reported rates of areal change for TNIC of up to −0.66 km 2 a −1 (197 to 169 km 2 ) between 1958 and 2000, while our own results give a nearly identical figure of −0.59 ± 0.03 km 2 a −1 (196.2 ± 9.9 km 2 to 173.2 ± 8.5 km 2 ) over this similar period.For GRIC, however, the shrinkage rate of −0.36 km 2 a −1 (135 to 120 km 2 ) reported by Sharp et al. (2014) over the period 1958-2000 is nearly 4 times more negative than our own figure of −0.10 ± 0.01 km 2 a −1 for 1952-1999 (131.8 ± 6.6 km 2 to 126.9 ± 6.3 km 2 ).Way (2015) recently reported that between 1973-1975 and 2010-2013, the area of TNIC decreased by 22 % (199.1 km 2 to 154.8 ± 7.5 km 2 ), while that of GRIC reduced by 18 % (134.3 km 2 to 110 ± 0.9 km 2 ).Although results slightly differ between these studies, our results agree within reported errors (when given).We hypothesize that those small disparities could be explained by the errors of interpretation due to snow patches around the ice caps, and/or by the different spatial resolutions and acquisition dates of the data used in the different studies (Paul et al., 2013).A comparison of the areal declines of GRIC and TNIC with those of other Baffin Island ice caps was already conducted in Way (2015) and is thus not presented here.
C. Papasodoro et al.: Changes of the two southernmost ice caps of the Canadian Arctic Archipelago Gardner et al. (2012) estimated that the average mass loss rate of all glaciers and ice caps on southern Baffin Island (South of 68.6 • N, excluding Penny Ice Cap) increased from −0.20 ± 0.05 m a −1 w.e. to −0.76 ± 0.12 m a −1 w.e.(i.e. a factor of 4) between the periods 1957-2006 and 2003-2009.This acceleration is more than twice that estimated over similar periods for northern Baffin Island glaciers (North of 68.6 • N, excluding Barnes Ice Cap).Barnes Ice Cap itself, located on central Baffin Island at elevations between 400 and 1100 m a.s.l., recently experienced a strong thinning acceleration (Sneed et al., 2008;Dupont et al., 2012), resulting in a mass loss rate of −1.08 ± 0.12 m a −1 w.e. between 2005 and 2011 (Gardner et al., 2012).The estimated mass loss rate on Penny Ice Cap between 2003 and 2009 is lower, at −0.52 ± 0.12 m a −1 w.e., a difference which Gardner et al. (2012) attribute to its much higher elevation range (up to ∼ 2000 m a.s.l.).The comparatively large mass loss rates experienced by GRIC and TNIC in the past half-century can be ascribed to differences in size and to the hypsometry of the ice caps, but also possibly to local climatic factors, as described below.

Regional context and climatic factors
The accelerating recession of glaciers and ice caps across the CAA in recent decades, and the concurrent increase in surface melt over the Greenland Ice Sheet, have been ascribed to a sustained atmospheric pressure and circulation pattern that favors the advection of warm air from the northwest Atlantic into the eastern Arctic and over western Greenland (Sharp et al., 2011;Fettweis et al., 2013).This situation has led to warmer, longer summer melt periods on glaciers of the eastern CAA, and this largely accounts for their increasingly negative mass balance (Weaver, 1975;Hooke et al., 1987;Koerner, 2005;Sneed et al., 2008;Gardner and Sharp, 2007;Gardner et al., 2012).
In the southern Baffin Island region, annual and seasonal mean air temperatures have generally increased since 1948 (except in the spring), but not monotonically (Vincent et al., 2015).At Iqaluit, seasonal trends from 1948 to ∼ 1990 were non-significant or slightly negative (spring).Thereafter, temperatures rose, particularly in autumn (SON; +0.8 • C decade −1 ) and winter (DJF; +2.9 • C decade −1 ), both of these trends being significant at the 95 % level (p < 0.05), even when autocorrelation is accounted for (Ebisuzaki, 1997).Climate records from stations further south (e.g.Resolution Island, 61.5 • N) are unfortunately too discontinuous to allow quantification of temperature trends on Meta Incognita Peninsula, but these are probably close to those observed in Iqaluit.Although GRIC and TNIC are only separated by 17 km, they did not experience the same historical and recent rates of shrinkage and mass loss (see the Sect.5), and part of the difference is likely due to differences in hypsometry, which strongly influences the response of glaciers to a given climate forcing (Oerlemans et al., 1998;Davies et al., 2012;Hannesdóttir et al., 2015).GRIC lies at a slightly higher altitude than TNIC, with 77 % of its area above 600 m a.s.l., compared to 68 % for the TNIC, and is therefore expected to have a slightly less negative mass balance, as our observations confirm.
A factor that may have indirectly contributed to the accelerating rate of glacier recession on southernmost Baffin Island is the decline in summer sea ice cover in this region (Fig. 9b), one of the steepest observed across the entire CAA (up to −16 % decade −1 since 1968; Tivy et al., 2011).In the Hudson Bay, Hudson Strait, Foxe Basin region, up to 70-80 % of the sea ice decline since 1980 has been attributed to warmer spring and/or autumn surface air temperature, wind forcing accounting for the balance (Hocheim and Barber, 2014).The retreating sea ice cover in Hudson Strait, immediately south of Meta Incognita Peninsula, has been accompanied by a particularly large rise in surface air temperature during autumn months (SON) during or after the sea ice cover minimum, and the rate of autumn warming between 1980 and 2010 (0.15 • C a −1 ) is estimated to have been three times greater than the mean between 1950 and 2010 (0.05 • C a −1 ; Hocheim and Barber, 2014).A consequence of the sea ice retreat in this sector has been to increase the net solar flux absorbed annually at the sea surface at an estimated average rate of ≥ 0.8 W m −2 a −1 over the period 1984-2006, and probably faster after the mid-1990s when sea ice decline accelerated (Matsoukas et al., 2010;Hocheim and Barber, 2014).Meanwhile, the temperature record from Iqaluit (Fig. 9a) indicates that while the cumulative sum of positive degree-days (PDD) during the glacier ablation season (April-November) was relatively constant prior to 1990 (no significant trend), it increased at a rate of ∼ 3.8 degreeday a −1 (p < 0.05) after 1990.The clearest increases in PDD occurred between summer and autumn (June to October: 0.24 to 1.46 degree-day a −1 ; p < 0.05), while trends in spring months (April and May) were comparatively very slight or negligible.These observations suggest that while rising summertime temperature undoubtedly remain the main driver for the mass losses of GRIC and TNIC over recent decades (Gardner et al., 2012), the annual mass loss rate could be enhanced by a lengthening of the melt season into the autumn linked to the later freeze-up in Hudson Strait (Hocheim and Barber, 2014).The early autumn weeks, in particular, are a period during which ice-free, open waters surrounding Meta Incognita Peninsula are a large local source of heat to the lower troposphere, while the frequent low-level cloud cover in this season may contribute a further downwelling longwave radiation flux (Dunlap et al., 2007).In this respect, the situation for GRIC and TNIC could differ from that of Barnes Ice Cap (70 • N) on central Baffin Island, where the lengthening of the melt season has been attributed to more frequent early spring thaw events (Dupont et al., 2012).Spaceborne monitoring of the melt duration over GRIC and TNIC by passive microwave sensing would help to verify if these two glacierized sectors of Baffin Island respond to regional warming in different ways.

Conclusions
This paper highlighted historical and recent trends in area, elevation and mass changes for the two southernmost ice caps of the Canadian Arctic Archipelago, Grinnell and Terra Nivea Ice Caps.Our analysis is based on multiple data sets and uses an original approach where ground control points for the photogrammetric processing of old aerial photographs are derived from sub-meter resolution Pléiades satellite stereo-images.This approach takes full advantage of the highly precise Pléiades products and represents an important advance for eventually unlocking the vast archives of historical aerial photographs.
Results show that the areal extent of TNIC is 34 % smaller in 2014 when compared to the end of the 1950's extent, while GRIC shrank by nearly 20 % between 1952 and 2014.Both ice caps also experienced an acceleration of their shrinkage rates since the beginning of the 21st century.
The historical glacier-wide mass balance for GRIC was estimated to be −0.37 ± 0.21 m a −1 w.e.  and slightly more negative for TNIC at −0.47 ± 0.16 m a −1 w.e.(1958/59-2014).Between 2007 and 2014, the mass balance of TNIC was of −1.77 ± 0.36 m a −1 w.e., a rate 5.9 as negative as the mass balance of −0.30 ± 0.19 m a −1 w.e.measured between 1958/59 and 2007.This is also twice as negative as the average mass balance obtained between 2003 and 2009 for other larger ice caps in the southern part of Baffin Island (Gardner et al., 2012).
The 2007-2014 mass balance of TNIC is among the most negative multi-annual glacier-wide mass balances measured to date, comparable to other negative values observed in the southern mid-latitudes (e.g.Willis et al., 2012;Berthier et al., 2009) or in southeast Alaska (Trüssel et al., 2013).Given the absence of calving glaciers for TNIC, its high rate of mass loss can only be explained by negative surface mass balance due to an ELA that, for most years, is above the maximum ice cap altitude.Nonetheless, this similarity in rate of mass loss underlines the strong sensitivity of maritime low-elevation ice bodies to the currently observed climate change at midlatitudes and in polar regions (Hock et al., 2009).The recent acceleration of ice cap wastage on Meta Incognita Peninsula is linked to a strong near-surface regional warming and a lengthening of the melt season into the autumn that may be reinforced by sea ice cover reduction and later freeze-up in Hudson Strait and nearby marine areas.

Figure 2 .
Figure 2. (a) Pléiades orthoimage of Grinnell Ice Cap (3 August 2014) superimposed with areal extents from 1952, 1999 and 2014.(b) Pléiades orthoimages of Terra Nivea Ice Cap (14 August 2014 on the east side and 26 August 2014 on the west side) superimposed with areal extents from 1958/59, 1998, 2007 and 2014.The overlapping area between the two orthoimages is represented by the black dashed polygon.(c) Historical and recent area changes for both ice caps.Error margins were estimated at 5 % for historical areas and at 3 % for 1998 and later outlines (Sect.4.4.3.).

Figure 3 .
Figure 3. Elevation change rates (m a −1 ) on Grinnell Ice Cap between August 1952 (Aerial photo DEM) and August 2014 (Pléiades DEM).For this figure as well as for the next ones, no color (i.e.hillshade is visible) represents no data.

Figure 4 .
Figure 4. Elevation change rates (m a −1 ) on Grinnell Ice Cap between March/April 2004 (ICESat and in situ GPS points) and August 2014 (Pléiades DEM).Bottom right graph shows historical (1952-2004) and recent (2004-2014) rates of elevation changes along the 203 points contiguous with the 1952 DEM (represented as black dots on the map).

Figure 8 .
Figure 8. Representation of the same geomorphological feature on ice-free terrain surrounding GRIC using three different technologies, namely an aerial photography (August 1952), a Pléiades panchromatic band (3 August 2014) and a Landsat 8 panchromatic band (15 August 2014).Note the very fine resolution of the Pléiades panchromatic band (70 cm), in comparison to the Landsat 8 panchromatic band (15 m), allowing to retrieve bedrocks and ice-free features on archives aerial photos and thus to collect GCPs (e.g. at the bedrock localised by the yellow arrow).

Figure 9 .
Figure 9. (a) Annual anomalies in total positive degree-days (PDD) recorded from April to November at the Iqaluit weather station, 1952 to 2014, based on Homogenized Canadian Historical Climate Data (Vincent et al., 2015).(b) Anomalies in total sea-ice covered area during the summer and autumn (25 June-19 November) over Hudson Strait (full line) and Davis Strait (dashed line), 1968-2014.Data provided by the Canadian Ice Service.For region boundaries, see Tivy et al. (2011), Fig. 4.

Table 2 .
Historical and recent glacier-wide mass balances for both ice caps.