Remote Sensing of Environment Rigorous 3D change determination in Antarctic Peninsula glaciers from stereo WorldView-2 and archival aerial imagery

This paper presents detailed elevation and volume analysis of 16 individual glaciers, grouped at four locations, spread across the Antarctic Peninsula (AP). The study makes use of newly available WorldView-2 satellite stereo imagery to exploit the previously untapped value of archival stereo aerial photography. High resolution pho- togrammetric digital elevation models (DEMs) are derived to determine three-dimensional glacier change over an unprecedented time span of six decades with an unparalleled mean areal coverage of 82% per glacier. The use of an in-house robust surface matching algorithm ensured rigorous alignment of the DEMs to overcome inherent problems associated with processing archival photography, most notably the identi ﬁ cation and correction of scale error in some datasets. The analysis provides insight into one of the most challenging and data-scarce areas on the planet by expanding the spatial extent north of the AP to include previously un-studied glaciers located in the South Shetland Islands. 81% of glaciers studied showed considerable loss of volume over the period of record. The mean annual mass loss for all glaciers yielded 0.24 ± 0.08 m.w.e. per year, with a maximum mass loss of up to 62 m.w.e. and frontal retreat exceeding 2.2 km for Stadium Glacier, located furthest north on Elephant Island. Observed volumetric loss was broadly, though not always, correlated with frontal retreat. The combined mass balance of all 16 glaciers yielded − 1.862 ± 0.006 Gt, which corresponds to − 0.005 mm sea level equivalent (SLE) over the 57 year observation period.


Introduction
The Antarctic Peninsula (AP) is one of the most inaccessible regions on Earth, with its remoteness and extreme weather conditions limiting human understanding of this fragile mountain glacier system. The AP is a complex system comprised of > 1590 glaciers ) that drain a narrow and high mountain plateau. The average mountain height is 1500 m, with the highest peaks rising to > 3000 m above sea level. As such, the region can be considered as a glaciated mountain range that differs fundamentally from the ice sheets covering East and West Antarctica. Most of the glaciers on the west and north-east coasts of the AP are marine-terminating, either as tidewater glaciers or with a short floating portion. The AP occupies < 1% of the entire grounded Antarctic ice sheet, but has the potential to significantly contribute to sea-level change (Meier et al., 2007;Pritchard and Vaughan, 2007;Radic and Hock, 2011). The AP ice sheet itself accounts for 25% of all ice mass losses from the Antarctic region, despite comprising only 4% of the continental area (Shepherd et al., 2012). Simulations of future sea level change in the AP suggest that omission of tidewater glaciers could lead to a substantial underestimation of the ice-sheet's contribution to regional sea level rise (Schannwell et al., 2016).
Similar to other large mountain near-polar glacier systems such as Alaska and Patagonia (Barrand et al., 2013a), the AP is most likely influenced by both atmospheric and oceanic changes (Cook et al., 2016;Vaughan et al., 2001). Its sensitivity to warming has manifested in many ways, one of the signs being the progressive retreat and subsequent collapse of numerous ice shelves such as Larsen A in January 1995, Wilkins in March 1998, and Larsen B in February/March 2002(Cook and Vaughan, 2010Scambos et al., 2004). In total, seven out of 12 ice shelves around the AP have either almost completely disappeared or retreated significantly in the second half of the 20th Century (Cook and Vaughan, 2010;Fox and Vaughan, 2005). Increased ice T velocities (Pritchard and Vaughan, 2007;Scambos et al., 2004;Seehaus et al., 2015;Wang et al., 2016), as well as the retreat and thinning of glaciers and ice caps (Cook et al., 2005;Ferrigno et al., 2006;Scambos et al., 2014) have also been reported.
Atmospheric warming of the AP in the second half of the 20th Century is well documented. Vaughan (2006), for example, reported a 74% increase in the number of positive degree days at the Faraday/ Vernadsky Station (65°15′ S, 64°16′ W) between 1950 and 2000. According to Sobota et al. (2015) the air temperature on King George Island increased by 1.2°C between 1948 and 2012. The mean temperature trend for AP stations between 1949 and 2002 was reported by Jacka et al. (2004) to be approximately + 4.4°C per century. However, according to Turner et al. (2016), this warming process has reversed in the 21st Century. Turner et al. (2016) identified mid-1998 to early-1999 as the turning point between warming (+ 0.32 ± 0.20°C per decade, 1979-1997) and cooling periods (−0.47 ± 0.25°C per decade, 1999-2014) and attributed the decadal temperature changes to the extreme natural internal variability of the regional atmospheric circulation, rather than to the drivers of global temperature change. This does not necessarily mean that the retreat of glaciers will slow down or stop. Another recent study by Cook et al. (2016) argues that the primary cause of glacier retreat in the AP is ocean-led rather than atmospheric-driven.
Due to a lack of detailed observations, many global mountain glacier mass balance inventories either do not include the AP (Dyurgerov, 2002) or use proxy values, such as the global average, instead (Leclercq et al., 2011). A recent glacier basin inventory by Cook et al. (2014) provides information on 1590 AP glacier basins in the form of frontal and areal changes. The study showed a north-south gradient of increasing ice loss across the AP, with 90% of 860 marine-terminating glaciers shown to have reduced in area. This study was, however, compiled based on a relatively coarse 100 m ASTER DEM. A recent surface mass balance study of the AP over the period 1979-2014 by van Wessem et al. (2016) estimated mass change at 5.5 km resolution based on the regional atmospheric climate model RACMO2.3 and a firn densification model (Ligtenberg et al., 2011). The average AP icesheetintegrated surface mass balance, including ice shelves, was estimated at 351 Gigatonnes/year with an inter-annual variability of 58 Gigatonnes/ year, and the western AP dominating the eastern AP. However, to this day, there remains a lack of detailed and high resolution measurements relating to almost all AP glaciers, and in particular there is limited information on mass balance change. Radic and Hock (2011) estimated the total sea-level rise from global mountain glaciers by 2100 as 0.124 m ± 0.037 m, with the most significant contribution from glaciers in Arctic Canada, Alaska and Antarctica (excluding ice sheets). Schannwell et al. (2016) predicted a contribution to sea level rise from AP tide water glaciers and ice-shelf tributary glaciers (split equally) in the order of 0.028-0.032 ± 0.016 m by 2300. Pritchard and Vaughan (2007) assessed that the annual sea level contribution from the AP region increased by 0.047 ± 0.011 mm between 1993 and 2003. Shepherd et al. (2012) reported ice sheet mass balance of the AP between 1992 and 2011 as − 20 ± 14 Gigatonnes/year. Hence, while the AP is believed to be a considerable component of the overall Antarctic ice imbalance (Shepherd et al., 2012), denser spatial sampling is required to better understand its contribution to sea-level rise.
The validation of predictive ice-loss models, and the consequent potential contribution to sea-level, requires accurate understanding of historical change and its drivers. It is therefore crucial to accurately estimate mass loss, to identify temporal and spatial patterns, and establish whether changes may be counter-balanced by increased snowfall and snow accumulation (Kunz et al., 2012a;Nield et al., 2012). Due to the rough terrain and inaccessibility to ground-based survey, the only practical source of information that can provide data over extended regions of the AP with relatively high temporal resolution is remote sensing. However, despite the increasing ubiquity of satellite observations, spatio-temporal records are generally too sparse and recent to enable the identification of long-term trends and variations (Fox and Cziferszky, 2008).
This study improves existing records by presenting new information for sixteen glaciers distributed across the western coast of the AP. This is carried out by utilising modern-day stereo satellite images and aerial photographs dating back to the 1950s from largely unexplored archives of > 30,000 frames held by British Antarctic Survey (BAS) and the United States Geological Survey (USGS) (Fox and Cziferszky, 2008). The study builds on existing literature that has attempted to quantify changes in individual AP glaciers (Cook et al., 2005;Fox and Cziferszky, 2008;Kunz et al., 2012a;Kunz et al., 2012b;Seehaus et al., 2015). These previous studies have focussed largely on observations made at the glacier fronts (Cook et al., 2005;Kunz et al., 2012a), primarily as a result of limited coverage of archival photography as well as poor image texture and contrast, or were time-limited to the last two decades (Seehaus et al., 2015). This research quantifies elevation changes across more complete glacier extents, at higher resolution than before (submetre image pixels of WorldView-2 and aerial photography) and over a time span of almost six decades. Furthermore, the spatial extent of the studied changes has been expanded to the north of the AP to include glaciers located in the South Shetland Islands.
In contrast to Kunz et al. (2012a) and Miller et al. (2009), both of which used Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) 30 m gridded DEMs, this study utilised sub-metre resolution WorldView-2 imagery (DigitalGlobe, 2016). In the longterm, sub-metre spatial resolution satellite imagery such as that provided by WorldView-2 has the potential to independently shift understanding of the processes taking place in the AP. Although at this moment temporal records are too recent to allow such imagery to be the only source for reliable long-term change studies, the data can facilitate the accurate registration of datasets from alternative sources. A recent study by Wang et al. (2016) presented the use of WorldView-2 for registration and ortho-rectification of declassified ARGON images. WorldView-2 images have been successfully used in several studies of glacier change around the world (Chand and Sharma, 2015;Karimi et al., 2012;Osipov and Osipova, 2015;Racoviteanu et al., 2015;Wang et al., 2016;Yavaşlı et al., 2015), however, in most cases only planimetric observations were used in analysis rather than rigorous 3D photogrammetric observations from stereo-imagery, as presented here. A full photogrammetric workflow allows for a high degree of quality control in the process of image orientation, as well as provides the opportunity for 3D validation of results e. g. manual verification and editing of DEMs in a stereovision environment. Such opportunity is eliminated in the case of 2D ortho-imagery as the third dimension is suppressed.
This study follows the proof-of-concept research presented in a pilot study by Fieber et al. (2016), where surface matching techniques were used to demonstrate elevation changes in three glaciers around Lindblad Cove, AP. DEMs were generated from 2014 stereo WorldView-2 data and a corresponding block of 1957 aerial archival imagery. Here, further results following a refined methodology are presented for 16 glaciers and additional problems encountered in the processing of archival imagery, and their solutions, are discussed.

Study areas
Sixteen glaciers selected for this study are grouped in four sites: Elephant Island (EI), King George Island (KGI), Lindblad Cove (LC) and Anvers Island (AI). The glaciers are spread across latitudes of 61°to 64°S outh and longitude of 54°to 63°West. The fact that the glaciers are clustered at four locales limits their spatial extent compared with 16 glaciers at different locations, but at the same time allows observation of whether glaciers of similar size in very similar topographic settings exhibit the same behaviour. Fig. 1 shows the location of the study sites  (Bindschadler et al., 2008;LIMA, 2015).  Fieber et al. (2016), these datasets were subsequently reprocessed using refined methodology and are included here for completeness.

Archival photography
The source of archival imagery was the Falkland Islands and Dependencies Aerial Survey Expedition (FIDASE) conducted in austral summers 1955-56 and 1956-57. The photographs from these campaigns were obtained at nadir specifically for topographic mapping purposes and in the context of archival imagery of this era they are of relatively high quality with systematic regional coverage and good stereo overlap ( EROS is the remote sensing data management, systems development and research field centre of the USGS that provides digital scans of archival imagery. Most aerial imagery is digitised using a Phoenix V System at 25 μm resolution. The Phoenix V System uses a scanning back camera rather than a flatbed scanner and requires a 1.5 to 2-minute long exposure per frame to create a product similar to that obtained from a flatbed scanner (Smith and Longhenry, 2008). On enquiry, however, it transpired that the FIDASE images were scanned on a flatbed photogrammetric scanner (Zeiss or Leica depending on the dataset) due to the film suffering from vinegar syndrome (age related degradation). The images were scanned at 21 μm resolution and subsequently resampled to 25 μm to ensure consistency with other EROS datasets (private correspondence with R. Longhenry and T. Smith, September 2016). The ground sample distance (GSD) of FIDASE imagery scanned at 25 μm is approximately 0.6 m.

Satellite imagery
WorldView-2 Ortho Ready 2A Stereo Imagery products (DigitalGlobe, 2016) were purchased from DigitalGlobe covering the 16 glaciers in four separate scenes (four image stereo-pairs). WorldView-2 images are geo-referenced with demonstrated geolocation accuracy without ground control of < 3.5 m circular error, 90% confidence (LE90) (DigitalGlobe, 2016). The geocoding information is provided in the form of Rational Polynomial Coefficients (RPCs). Toutin et al. (2012) showed that the absolute elevation extraction accuracy from stereo WorldView-2 data without ground control points was 3.6 m (LE90) for bare terrain. Each stereo scene was acquired on a single pass of the WorldView-2 satellite over the AP. The image quality varied across the four scenes: two of the scenes (EI, AI) were overexposed, lacked texture and exhibited striping, while the remaining two were excellent with clear detail and good texture. The specific acquisition dates, are listed in Table 2. The nominal resolution of WorldView-2 images is 0.46 m at nadir for panchromatic images and 0.52 m at 20 o off-nadir. The specific maximum GSD for each of the four acquired stereo pairs is listed in Table 2. Table 2 shows that there is a substantial in-year time gap between the archival aerial photography and the recent satellite imagery for King George Island (77 days) and Lindblad Cove (80 days). The archive of stereo-satellite imagery was carefully searched for cloud-free images close to the corresponding Julian day of the archival photographs, but at the time of image selection stereo-coverage of Antarctica was far from complete. The acquisition of new satellite imagery was considered, but there was low confidence in the likelihood of attaining cloud-free imagery within the time-scale of the work. The northern AP is a consistently cloudy region with, for example, mean cloud cover of 7 oktas and 271 precipitation days per year (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) at Teniente Rodolfo Marsh Martin Station, King George Island, and mean of 6.5 oktas and 250 precipitation days  at Vernadsky Station, Argentine Islands (Kirchgaessner, 2010).

Seasonality of image sources
The substantive in-year time difference has the potential impact of skewing the change measurements due to different surface snow conditions. There are few data on annual snowfall in the northern AP, but Goodwin et al. (2016) describe an ice core from Bruce Plateau (66°S 64°W) that shows accumulation of 1.84 m/yr water equivalent (years 1900-2009). This is equivalent to about 6 m of snowfall in a year. Similarly, interpolated/modelled accumulation maps published by Turner et al. (2002) show 1.5-2.0 m and van Lipzig et al. (2004) 2.0-2.5 m.w.e. for the Lindblad Cove area.
The melt-season for the northern AP around Lindblad Cove starts in early to mid-October and continues through the austral summer until about mid-February (Barrand et al., 2013b), and is even longer on King George Island. Lindblad Cove satellite images were acquired on 4 April 2014, which is into the austral autumn, but examination of the images showed surface detail indicating that new winter snowfall had not yet started to accumulate. It is, however, possible that the mid-summer archival images from late December to mid-January hold approximately 2-3 m.w.e. of surface snow remaining from the winter after a period of melting from October onwards, coupled with compaction over time.

Methodology
Processing of archival imagery is often anything but a trivial task. Even when calibration information exists, which is often not the case, archived film may be damaged or degraded due to inappropriate storage and/or the passage of time and the radiometric quality of the sensors used for acquisition may not be very high (Fox and Cziferszky, 2008). Particularly in the case of areas covered by snow and ice, photography is often over-exposed due to extreme contrasts, resulting in a loss of surface texture. Poor viewing geometry and cloud cover can further limit the usability of archival photography. A rigorous frame selection process therefore needs to be applied and care needs to be exercised in the digitisation process if historic photography is to be used in metric studies. Glaciers for this study were selected on the basis of the quality of available archival data (e.g. minimal cloud cover, good texture) to ensure near-complete glacial coverage. These criteria are together a strong filter on the number of photographs suitable for processing (Table 3). For the selected glaciers with historical photograph coverage, a WV2 image of the best available quality was subsequently sought. Unfortunately, WV2 images suffer from similar constraints in terms of overexposure and cloud coverage and as such the quality of the imagery was not always ideal.
As opposed to recent developments of image-based fully-automated structure-from-motion measurement techniques, this study follows a rigorous photogrammetric approach. Structure-from-motion has its roots in the computer vision community and has recently gained a lot of interest in many geo-applications, including the cryosphere field (Nolan et al., 2015;Piermattei et al., 2015), for its ease of use and no need for specialist's knowledge. Nevertheless, despite the considerable benefits, such as no requirement for prior camera parameters knowledge, it is still prone to geometric errors and it often offers a 'black-box' approach which allows user's little or no control over the process. In the absence of ground control, coupled with the availability of camera calibration and classical stereo network geometry, the 'classical' photogrammetric approach was favoured in this study for its rigour and accuracy aspects. Both archival photography and modern satellite imagery was processed in Socet GXP software (BAE Systems, 2015) following a rigorous photogrammetric workflow that is based on the collinearity condition (see, for example, Kraus, 2007). The collinearity condition states that at the time of image acquisition, the position of the ground points in object space, their depiction in the image and the location of the sensor station all lie in a straight line. Each image is then considered as a bundle of rays that can be manipulated in a bundle adjustment. Here all bundles are translated, rotated and scaled in a single mathematical adjustment until corresponding rays from different bundles intersect at common ground points.
The workflow is illustrated in Fig. 2. At first an approach akin to that of Fox and Cziferszky (2008) and Barrand et al. (2009), in which pseudo-GCPs are extracted from the modern geocoded dataset and introduced into the bundle adjustment of the archival data, was adopted to ensure that the datasets were aligned within the pull-in range of the subsequent least squares processing procedure. Secondly, to minimise misalignment errors, a robust least-squares surface matching technique was applied. The in-house surface matching algorithm utilised here was originally developed by Mills et al. (2005) for coastal change monitoring and further enhanced by Miller et al. (2007) and (2008) for wider geohazard application. In glacier monitoring applications, various incarnations of the approach have been demonstrated by Miller et al. (2009), Kunz et al. (2012a) and Fieber et al. (2016).
Four WorldView-2 stereo-pairs, provided with orientation and geocoding information, were visually examined via stereo-viewing and if Y-parallax was observed, automatic point matching and relative orientation using full-covariance strategy was performed (resulting in RMSE < 0.2 pixel). GCP candidates were subsequently chosen via stereo-viewing of mountainous areas that were presumed to be stable over the period of observation.
Interior orientation of FIDASE images (see Table 3 for the number of selected frames for each site) was carried out by manual measurement of fiducial marks using available camera calibration information. It should be noted that interior orientation will largely deal with any systematic shrinkage or linear deformation of the film that has occurred over time (and which would therefore be present in the scanned digital imagery). The additional subsequent use of the 9-parameter least squares surface matching process would also help mitigate any residual error from such issues. Subsequent relative orientation was initially performed using the automatic tie point measurement strategy embedded in the software and refined manually with the help of stereoviewing. Finally, absolute orientation using full-covariance bundle adjustment strategy was carried out by the introduction of pseudo-GCP candidates transferred from WorldView-2 satellite data and their manual measurement in stereo (see Fig. 3 and Table 3 for the number of pseudo-GCPs at each site). Different acquisition angles, as well as different snow cover conditions coupled with glacier elevation change across six decades, rendered the process of selection and transfer of the pseudo-GCPs between two datasets time-consuming and non-trivial. Extreme care was taken throughout the process to ensure correct identification of points, even distribution across the site and to facilitate good initial alignment of the datasets.
Once the orientation of FIDASE images was deemed satisfactory, photogrammetric DEMs were extracted for both satellite and archival datasets using Socet GXP Next Generation Automatic Terrain Extraction (NGATE) low contrast strategy as TIN DEMs with approximately 1 m post spacing. Manual editing of the DEMs was subsequently performed in stereo-viewing to ensure removal of gross image-matching errors and to improve the coverage in the areas of little texture, where the automatic algorithm failed to produce satisfactory results. Rocky areas, considered to be stable over time, were manually delineated using WorldView-2 scenes and refined by comparison to archival data to avoid inclusion of areas of significant glacier change, where rock outcrops have been exposed over the period of time and were not visible in the early imagery. The surface matching algorithm of Kunz et al. (2012b) was then applied to 'stable area' DEMs to refine the alignment of the archival dataset to its modern day reference.
The least-squares algorithm iteratively minimises vertical (as used here) or Euclidean differences between points in the clipped stable areas on the floating surface (archival dataset) and corresponding areas on the fixed reference surface (modern data). A weighting function, based on a maximum likelihood estimator, is included in the algorithm to down-weight potential outlier points, such as those related to small areas of temporal change or simply gross errors, to avoid biasing the solution into an erroneous alignment (Miller et al., 2008). The software provides the choice of the number of transformation parameters to be retrieved. In this study, instead of a usual Helmert seven-parameter transformation (three translation, three rotation, and one scale parameters), a nine parameter transformation was enabled that incorporated different scales in each Cartesian axis.
Once the transformation parameters were successfully retrieved over stable terrain, the transformation was applied to the full archival DEM. The differences in elevation and volumetric changes between different epochs were then computed using LSS software (McCarthy, 2015). Further, mean elevation change as metre water equivalent (m.w.e.) over the period of record, as well as annually, was calculated based on the volume to area relationship (Fieber et al., 2016;Hock,  2010). Assuming that most of the lowering occurs at the density of ice, an ice density correction factor of 0.917 was applied here to express the rates in water equivalent units, (Paterson, 1994;Vaughan et al., 1994) as well as to convert volume to mass. Sea level equivalent was computed assuming an ocean area of 362.5 × 10 6 km 2 (IPCC, 2013). Further details of the methodology can be found in Fieber et al. (2016).

Archival imagery processing
The results of the initial orientation of the FIDASE imagery at all four sites, carried out in Socet GXP, are summarised in Table 4. In the cases of Elephant Island and Lindblad Cove, two separate projects were set up to orientate the photography. At Lindblad Cove this was due to two different image acquisition dates (see Table 2), while in the case of Elephant Island this was due to a problematic orientation process resulting from the high percentage of the FIDASE frames filled with sea and little texture in the overlap area between strips, necessitating the separate processing of the two flight lines. Orientation statistics indicate that the results were best at Lindblad Cove, while Anvers Island and Elephant Island South produced the highest RMSEs. This was likely the combined result of several factors: (1) the archival photography suffered from a scaling error, seemingly introduced at the digitisation stage as a result of resampling from 21 μm to 25 μm pixel size, which subsequently led to misinterpretations in the processing software and poor orientation results; (2) both Elephant and Anvers Islands had relatively poor quality WorldView-2 imagery (striping, lack of texture) which impacted on identification of pseudo GCPs; (3) archival imagery very often was taken in such a way that high percentages of the frames were filled with sea (particularly Elephant Island South), rendering frame orientation problematic; (4) in the case of Anvers Island, terrain elevation changed from 0 to > 2500 m across only five FIDASE frames, relative to a flying height of about 4000 m, causing considerable scale change, difficulty in matching tie points between images and impacting orientation results.

DEM quality assessment
The assessment of absolute DEM accuracy was not possible due to a lack of existing known control points in the areas being surveyed. Instead, Fieber et al. (2016) estimated internal quality of the DEMs generated by SOCET GXP through comparison of manual and automatic measurements at Lindblad Cove, making the assumption that manual observations were of higher quality. Points on the glacier surface were manually measured in a semi-gridded pattern in stereo, both in WV2 and FIDASE imagery for Lindblad Cove projects. The measurements were performed in ten randomly selected areas of interest, each approximately 250 by 250 m in size. Those measurements were subsequently compared to corresponding areas in WV2 and FIDASE DEMs generated automatically through image matching by Socet GXP. The standard deviation of the elevation differences was adopted as a quality  measure for WorldView-2 ( ± 1.82 m) and FIDASE imagery ( ± 2.57 m). The uncertainty in elevation differences between two surfaces was then calculated based on error propagation and yielded ± 3.14 m. This error estimate is applicable to this study.

DEM matching
Surface matching was applied to all archival DEMs (stable areas only) to improve their alignment with the corresponding WorldView-2 DEMs. Initially, the software was set to search for the seven parameters of a Helmert transformation. The match achieved at Lindblad Cove was exceptionally good (Fieber et al., 2016), however manual inspection, including visual checks of point cloud cross-sections, indicated that in the case of other sites, particularly Elephant and Anvers Island DEMs, the fit was unsatisfactory (elevation of sea level varied significantly when inspected in stereo). Even where surfaces appeared to statistically correctly co-register, checks on sea-level were outside acceptable error margins. Surface matching was therefore repeated, allowing independent scale changes in all three Cartesian axes through a nine parameter transformation. The Lindblad Cove dataset was also reprocessed using a nine parameter transformation to ensure consistency. Table 5 lists the determined transformation parameters, while Table 6 and Fig. 4 show elevation difference statistics and histograms, respectively, between stable areas of modern and archival DEMs for all study areas before and after application of nine parameter surface matching. Allowing differential scale changes in each Cartesian axis significantly improved the fit of the surfaces. The improvement in the elevation difference statistics post-match can clearly be seen in Table 5 for all sites. Visual comparison of cross-sections was also satisfactory with sea-level mostly within the uncertainty error of ± 3.14 m and outliers not exceeding three sigma. Surface matching can therefore be deemed successful in improving the alignment of the surfaces. Fig. 5 illustrates elevation differences for all sixteen glaciers while corresponding histograms of elevation changes are presented in Fig. 6. It should be noted that the histograms correspond to the points of measurement and may be therefore slightly biased towards the areas of glacier with higher point density (lower point density areas correspond to the manual elevation measurements, where image matching failed due to poor texture), they are however very illustrative in terms of overall glacier change. The statistical information is complemented by Table 7 which summarizes the volumetric and mass change in those glaciers over a time period of 56 to 58 years. Elevation change maps as well as the statistics of 13 glaciers show predominance of loss, while three remaining glaciers Emerald Icefalls A at KGI (Fig. 5J), McNeile at LC (Fig. 5O) and Schoeling at LC (Fig. 5N) exhibit accumulation. The glaciers that showed the highest volume loss are Stadium at EI (Fig. 5A), Lange at KGI (Fig. 5F), Znosko at KGI (Fig. 5D) and Admiralen at KGI (Fig. 5E). Stadium is one of three glaciers in this study that is located furthest north, on Elephant Island, while Lange, Znosko and Admiralen are situated on King George Island, also further north than Lindblad Cove and Anvers Island. The highest loss in the furthest north glaciers is in agreement with Kunz et al. (2012a), who observed higher rates of lowering in the northern AP than in the south, although this study did not include glaciers as far north as on South Shetland Islands. Two of the three glaciers that exhibited predominance of accumulation over loss (McNeile and Schoeling) are located at Lindblad Cove. This in turn corresponds to the low mean ocean temperatures between 1945 and 2009 (< − 1°C) and mostly positive glacier areal changes in that region (east side of Bransfield Strait) reported in Cook et al. (2016). Similarly considerable loss of volume in Kraus at AI (Fig. 5P) coincides with positive mean ocean temperatures (> +1°C) and negative areal glacier changes on Anvers Island (Cook et al., 2016).

DEM matching
The post-matching RMSE of DEM differences over stable (generally rough mountainous) terrain reduced considerably in comparison to the pre-match RMSE for all sites with the exception of Lindblad Cove, where little improvement was noted. On average post-match RMSE yielded ± 18.32 m, which is comparable to values previously reported in literature (Kunz et al., 2012a;Kunz et al., 2012b). The elevation differences included in the RMSE can be attributed to rough and steep terrain used for matching. The DEM quality over significantly flatter glacier surface is expected to be much closer to the propagated uncertainty of ± 3.14 m (Section 4.1.2). However, due to the fact that surface matching also impacts on the resultant elevation changes, the true accuracy of those changes will likely lie somewhere between the two accuracy estimates ( ± 3.14 m and individual post-match RMSE). Therefore, for the purpose of uncertainty in the mean estimation, the higher of the two accuracy values, i.e. post-match RMSE at 95%  Table 6 Statistics of elevation differences over stable terrain for all studied sites between FIDASE and WorldView-2 DEMs before and after surface matching using 9-parameter transformation. confidence interval, together with the number of observations at each site was used. Uncertainty in the net volume was then computed based on uncertainty in the mean and the glacier surface area. At the same time, surface matching helped highlight a problem with the processing. Table 4 summarizes Z-scale factor (S z ) determined by the software for all sites. In three out of four sites the Z-scale is considerably smaller than unity. Surface matching showed that Elephant and Anvers Island DEMs needed 11% Z-scale correction, while King George Island around 6.5%. Only in the case of Lindblad Cove was the Z-scale close to unity. That also explains the best orientation and the least improvement in the post-match statistics of the Lindblad site (Table 5). Such Z-scale error is usually indicative of a focal length calibration error. However, on checking, all of projects were found to have implemented camera parameters as per the appropriate calibration certificates. Nevertheless, the investigations revealed that the flying heights after absolute orientation in the Z-scale affected projects were much higher than the nominal flying altitude of 13,500 ft. Since the focal length was kept fixed, it appeared that the adjustment was compensating for an 'incorrect' focal length by increasing the flying height, indicative of a problem in image interior orientation.
In order to eliminate the possibility of error in the imagery/film itself, Elephant Island South negatives were re-scanned on the BAS Zeiss/Intergraph Photoscan TD scanner at 21 μm resolution. This imagery was subsequently orientated using the same pseudo GCPs as utilised for USGS-scanned images. The DEM was generated and equivalent stable areas used in nine-parameter surface matching with the same WorldView-2 DEM. BAS-scanned images orientated significantly better, although not ideally, indicating that other factors such as quality of WorldView-2 imagery and a significant part of FIDASE frames occupied by sea, played their role. Relative orientation RMSE of BAS-scanned images yielded 0.92 pixel (rather than 1.48 pixel) while absolute orientation RMSEs were 4.38 m in X, 10.34 m in Y, and 2.24 m in Z (in comparison to Elephant Island South using USGS imagery in Table 4 -9.27 m in X, 15.49 m in Y and 3.56 m in Z). Furthermore, there was no increase in determined flying heights and surface matching indicated no Z-scale variation (S z = 0.9969 rather than 0.8904).
Discussion with USGS/EROS confirmed that USGS scanned the images using Zeiss (KGI, EI, AI) and Leica (LC) photogrammetric scanners and that post-scan, the original scan resolution of 21 μm was resampled to 25 μm to ensure consistency with other EROS digital archive products. Seemingly, in the case of the Zeiss scans, this resampling process somehow corrupted the file headers and subsequent interior orientation through erroneous pixel size, which led to the scaling issue in the image derived products. This example shows that the processing of archival imagery needs to be treated with extreme care and that working with such imagery is far from a trivial task. The problem identified above also proves the validity of using rigorous orientation procedures in the processing of such imagery.

Glacial change
In all but two cases (Znosko at KGI (Fig. 5D) and Emerald Icefalls A at KGI (Fig. 5J)) the volume loss/gain was correlated with glacier frontal retreat/advance, respectively. Stadium atEI (Fig. 5A) retreated most significantly, by a maximum of − 2.25 km and lost a net volume of −0.819 ± 0.0014 km 3 with m.w.e. of − 62.1 (− 1.07 m.w.e. p.a.). Similarly, Lange at KGI (Fig. 5F) retreated by 1.45 km at maximum, losing a volume of − 0.785 ± 0.0003 km 3 with m.w.e. of − 50.2 (− 0.89 m.w.e. p.a.). The retreat (at maximum) of the remaining eleven glaciers varied between −70 m and −450 m. Only three out of the sixteen glaciers studied advanced: Znosko at KGI (Fig. 5D) (c. +100 m at maximum), McNeile at LC (Fig. 5O) (c. +130 m at maximum), and Schoeling at LC (Fig. 5N) (c. + 200 m at maximum). In the two cases where the observed volume loss/gain was not correlated with a respective retreat/advance, Znosko (KGI) advanced locally by almost 100 m over the 56-year period, but at the same time showed considerable elevation lowering, with m.w.e. of − 31.6 (− 0.56 m.w.e. p.a.). Conversely, despite exhibiting predominance of accumulation with m.w.e. of + 1.8 (+0.03 m.w.e. p.a.), Emerald Icefalls A at KGI (Fig. 5J), retreated locally by a maximum of 70 m. Both glaciers were studied over either their full (Znosko) or almost full (Emerald Icefalls A at 84%) areal coverage, therefore it is unlikely that the outcome would change with more extended treatment and these examples thus show that sometimes small frontal change, be it advance or retreat, is not always indicative of associated glacier volume/mass change. Fig. 7 summarizes the overall change in all studied glaciers over the period of record expressed in the form of m.w.e. p.a., while Fig. 8  K.D. Fieber et al. Remote Sensing of Environment 205 (2018) 18-31 presents elevation change rates per year over elevations above sea level (second order polynomial fit) for all glaciers grouped according to their location. The measurements for the latter plot were obtained along each glacier centreline and subsequently combined for each study site. The general pattern of surface lowering and frontal retreat appears consistent. Elephant Island and King George Island show the highest lowering rates towards lowest elevations which is related to the plots being biased towards two glaciers with the highest volume loss: Stadium and Lange, respectively. Lindblad Cove is the only site which exhibits overall accumulation stability across all elevationsno lowering at low altitudes is observed. As noted in Section 2.2.3, the fact that the WorldView images were acquired 80 days later in the year than the FIDASE photographs may influence the observations. In three study areas, namely Elephant Island, King George Island and Anvers Island, the lowering decays to zero within elevations of 400 m to 500 m above sea-level ( Fig. 8). At higher altitudes accumulation or overall stability can be noted. This is consistent with the findings of Kunz et al. (2012a), who reported the lowering to be limited to the regions below 400 m. In a study of mass balance changes for the King George Island ice cap, Rückamp et al. (2011) Rückamp et al. (2011) study is much shorter, the 0.20 m p.a. at 270 m above sea level is consistent with this study's findings for King George Island (Fig. 8). At the same time, the observed 1.4 m p.a. elevation change at 40 m is a much higher estimate of lowering towards the front of the glacier than presented here, indicative of an accelerating surface lowering rate in recent decades. The mean total lowering for all 16 glaciers over an average period of 57 years was approximately 14.0 m ± 4.8 m.w.e. (standard error of the mean -SEM) with a maximum lowering of up to 62 m.w.e. for Stadium Glacier. The mean annual lowering for all glaciers yielded 0.25 ± 0.08 m.w.e. p.a. This is in agreement with observations reported in Kunz et al. (2012a) of a mean frontal lowering rate of 0.28 ± 0.03 m.w.e. p.a. over an average period of 37 years. Kunz et al., however, determined mean lowering using a 1 km 2 sample area of each studied glacier within 1 km of the front, whereas here the complete or near-complete surface area of each glacier was used in the calculation. The average spatial coverage of all glaciers obtained in this study was 82% ± 5%. This represents a relatively good achievement in comparison to Kunz et al. (2012a), where the focus was on the glacier fronts and where glacier coverage was on average 20%, up to a maximum of 33%. The spatial distribution of changes across the AP indicates overall greater lowering rates in the north than in the south, agreeing with similar findings by Kunz et al. (2012a). Nevertheless, there is clear local variability in the change rate between neighbouring glaciers which may be related to the terrain topography, type of glacier, atmospheric and oceanic warming and other non-trivial factors. The fact that surface change across studied glaciers does not present a clear spatial pattern suggests that the processes driving these changes in the AP are not simple and cannot be easily generalised.
Data in Fig. 8 and Table 7 reveal that Elephant Island and King George Island exhibit the greatest combined mass loss of 0.809 ± 0.0028 Gt and 0.904 ± 0.0016 Gt, respectively, over the period of 57 years; Anvers Island (where only one glacier was observed) mass balance is also negative: − 0.212 ± 0.0009 Gt, while Lindblad Cove was the only site to show mass gain (0.063 ± 0.0009 Gt) in the study period. The combined mass balance of all 16 glaciers yielded − 1.862 ± 0.006 Gt (− 0.033 ± 0.0001 Gt p.a.), which corresponds to − 0.005 mm sea level equivalent (SLE) over the 57 year observation period. This value is close to the lower boundaries of the annual sea level rise contribution from the entire AP, which was computed based on runoff rates in 2000 by Vaughan (2006). Based on ICESat observations, Gardner et al. (2013) reported an annual change estimate of − 6 ± 10 Gt in the Antarctic and Sub-Antarctic (not including the northern AP) for the period of 2003-2009. The combined results of the 16 glaciers presented here constitutes 0.5% of that estimate.
Finally, Fig. 9 illustrates mass budget for all four study sites, in units of kg/m 2 per annum, overlaid on the regional glacier mass budget for the Antarctic and Sub-Antarctic, as presented by IPCC (2013). Although, the IPCC's's (2013) mass budget does not include the northern part of the AP (Gardner et al., 2013), where most of the sites studied  Table 7 Volumetric change statistics for sixteen studied glaciers. Ice density value of 0.917 Gt/km 3 used for volume to mass conversion. Sea level equivalent (SLE) computed using an ocean area of 362.5 × 10 6 km 2 for conversion (IPCC, 2013 here were located, the site-level estimates cover the span of mass budget modelled with climate data , fitting well with the IPCC's findings. The combined mass budget yielded − 251.9 kg/ m 2 per annum.

Conclusions
This study has determined high resolution elevation change for sixteen individual glaciers distributed across the AP over a period of almost six decades. In comparison to previous studies (Kunz et al., 2012a), the emphasis here was on achieving high glacier areal coverage (up to 100%), high resolution of spatial sampling (nominal 1 m DEMs) and long temporal data span for surface elevation comparison (~57 years). Although, 100% coverage of glacier surface was not always possible to achieve, due to lack of coverage by one or other of the datasets or because of clouds obscuring the surface or poor image Fig. 7. Mean annual metre water equivalent for all studied glaciers. Observation period for all glaciers is approximately 57 years. Background image -LIMA mosaic (Bindschadler et al., 2008;LIMA, 2015). K.D. Fieber et al. Remote Sensing of Environment 205 (2018) 18-31 texture, the glacier surfaces were sampled on average at 82% ± 5%. The lowest coverage of 44% and 48% was for Emerald Icefalls C (KGI) and Lange (KGI), respectively. This study has highlighted the potential and pitfalls of using archival aerial stereo-photography combined with modern imagery such as WorldView-2 stereo products for determining glacier elevation change in remote regions such as the AP. The use of archival and modern imagery with analogous spatial resolution, advances in photogrammetric processing capabilities and the use of a rigorous orientation workflow has facilitated accurate comparison of archival and modern DEMs of glacier surfaces. The surface matching technique proved particularly valuable in improving the fit between archival and modern surfaces and therefore making the elevation and volumetric change analysis both feasible and reliable. Surface matching helped to identify problems with archival imagery, most likely induced by resampling of the original scans, and corrected them using a nine parameter transformation. This highlighted the need for care and caution when working with archival photography, as well as the value of extending the employed surface matching solution to include independent scale parameters. Most importantly however, this emphasised the validity of using rigorous approaches, such as least squares surface matching, for co-registration when undertaking precise change detection.
Thirteen out of the sixteen glaciers studied showed elevation lowering and predominance of volumetric loss over gain. Volumetric loss was largely correlated with frontal retreat, however in two cases the observed volume loss/gain was not correlated with a respective retreat/ advance, highlighting that sometimes observed frontal change is not always indicative of associated glacier volume/mass change. The greatest change was observed at Stadium Glacier (EI), located furthest north, where frontal retreat exceeded 2.2 km (at maximum) and the volume reduced by − 0.819 ± 0.0014 km 3 (66% of the glacier surface studied) with m.w.e. of −62 m over 58 years. The second most significant change was observed at Lange Glacier (KGI) which retreated by over 1.4 km and lost volume of −0.785 ± 0.0003 km 3 (48% of the glacier surface studied) with m.w.e. of − 50 m over 56 years. Three out of sixteen glaciers, namely Schoeling (LC), Emerald Icefalls A (KGI) and McNeile (LC), showed predominance of accumulation over loss. Two out of those glaciers are located on the mainland, on the west coast of the AP rather than on an island and their location corresponds to the areas with low mean ocean temperatures reported in Cook et al. (2016). Accumulation, however, was not as significant as the loss observed in other glaciers, with m.w.e. not exceeding +2.8 m.w.e. over 57 years. A further two glaciers, Helava (EI) and Landau (LC), showed small losses with m.w.e. not exceeding −1.5 m.
The combined mass balance of all sixteen glaciers yielded − 1.862 ± 0.006 Gt, which corresponds to − 0.005 mm sea level equivalent (SLE) over the observation period of 57 years. Mass budgets of glaciers grouped into four sites expressed in kg/m 2 p.a. fitted well into the findings reported in IPCC (2013), covering the span of climatemodelled data. The observed surface lowering for glaciers grouped at four sites was confined to elevations below 400-500 m ASL for three out of four locations (King George Island, Elephant Island and Anvers Island), while Lindblad Cove exhibited overall accumulation stability across all elevations. The observed rates of surface lowering were greater in the northern AP than in the south, however high local variability was also observed, suggesting that the processes driving these changes in the AP are not simple and cannot be easily generalised.