Rapid demise and committed loss of Bowman Glacier, northern Ellesmere Island, Arctic Canada

Abstract

Using historical and recent aerial photography and structure from motion (SfM) multiview stereo (MVS) techniques, we reconstruct the 1959 and 2018 ice surface topography and determine the geodetic mass balance of Bowman Glacier, a small mountain glacier on northern Ellesmere Island. This is combined with optical satellite imagery to reconstruct the evolution in extent of the glacier over six decades, and ground-penetrating radar measurements of ice thickness to estimate the remaining ice volume. Between 1959 and 2020, Bowman Glacier lost 78% of its extent (reducing from 2.75 to 0.61 km²), while average annual area loss rates have nearly tripled in the past two decades. Over the 1959–2018 period, glacier-wide ice-thickness change averaged −22.7 ± 4.7 m, corresponding to a mean specific annual mass balance of −347.0 ± 71.4 mm w.e. a⁻¹. Projecting rates of area and volume change into the future indicates that the glacier will likely entirely disappear between 2030 and 2060. This study demonstrates the potential of SfM-MVS processing to generate elevation products from 1950/60s historical aerial photographs, and to extend observations of ice elevation and glacier volume change for the Canadian Arctic, prior to the satellite record.

Introduction

Glaciers and ice caps in the Canadian Arctic Archipelago (CAA) cover ~146 000 km², the largest ice-covered area outside of the ice sheets. Ice loss in the region has accelerated in recent decades due to enhanced meltwater runoff (Gardner and others, 2011; Lenaerts and others, 2013; Noël and others, 2018), with losses for the 2000–19 period occurring at one of the highest rates in the world at −57.1 ± 9.1 Gt a⁻¹ (split between −30.6 ± 4.8 Gt a⁻¹ for CAA North and −26.5 ± 4.3 Gt a⁻¹ for CAA South; Hugonnet, 2021). Regional mass-balance estimates have generally been derived from satellite altimetry, gravimetry and stereoscopic imagery, or regional climate models (Gardner and others, 2011; Lenaerts and others, 2013; Noël and others, 2018; Ciracì and others, 2020; Hugonnet, 2021). While valuable for monitoring ice-volume change on a large scale, these products generally suffer from limited spatial resolution and struggle to resolve changes on small ice masses with complex geometries. Additionally, small glaciers are notoriously underrepresented in both regional and global glacier inventories (Parkes and Marzeion, 2018). Meanwhile, smaller glaciers in the CAA are experiencing greater retreat rates and overall area loss than large ice masses, with some small glaciers disappearing entirely since 2000 (Thomson and others, 2011; Sharp and others, 2014; Serreze and others, 2017; White and Copland, 2018). For example, White and Copland (2018) reported that glaciers <1 km² across northern Ellesmere Island lost area at an average rate of 30.4% decade⁻¹ over the period 1999–2015, compared to a mean loss rate of 0.4% decade⁻¹ for glaciers >1000 km².

Determination of past glacier fluctuations requires long-term monitoring programmes, but regular observations of changes in ice extent, surface elevation and ice thickness exist for only a few individual glaciers across the CAA. In situ measurements are time-consuming, logistically expensive and are understandably impossible to implement on more than a limited sample of glaciers. Yet, the current trend of accelerated mass loss highlights the importance of long-term datasets for detecting change over decadal timescales, and the need for elevation products for deriving geodetic mass-balance estimates prior to the satellite record. Air photo surveys conducted by the Royal Canadian Air Force in the late 1950s and early 1960s provided the basis for the regional topographic maps for the CAA, which were then used to generate the first version of a digital elevation model (DEM) for the region, the Canadian Digital Elevation Data (CDED). However, earlier elevation products generated using photogrammetric methods contained considerable errors, particularly in the accumulation area of glaciers and ice caps with few bedrock features to match. For example, on Axel Heiberg Island in the northern CAA, the CDED contained elevation errors of >90 m (Cogley and Jung-Rothenhäusler, 2004), making it essentially impossible to derive geodetic mass-balance estimates. To our knowledge, very few recent studies have made direct use of photogrammetry techniques from 1950/60s aerial photography to derive elevation data products over glaciers in the region (Burgess and Sharp, 2008).

More recently, the use of digital aerial photography and structure from motion (SfM) photogrammetry coupled with multiview stereo (MVS) reconstruction techniques has become common in glaciological research. In the Canadian Arctic, SfM-MVS processing has been applied to a handful of glaciers for mapping present-day surface topography (Thomson and Copland, 2016) and supraglacial drainage networks (Bash and others, 2022), and quantifying short-term surface ablation (Bash and others, 2018). Thomson and Copland (2016) applied SfM-MVS methods to oblique air photos from a 2014 survey, creating a DEM of the surface of White Glacier in Expedition Fiord, Axel Heiberg Island, and determined geodetic mass change from comparison with a 1:10 000 topographic map derived from 1960 aerial photography (Thomson and others, 2017). Contours in the 1960 map were digitised to produce a DEM of the area which was then compared to the 2014 DEM (Thomson and others, 2017). The downside of this method is that manual digitisation is time-consuming and impractical for large study areas, and errors in topographic map data and uncertainties associated with interpolating elevation contours are propagated to the DEM and to surface change measurements (James and others, 2012). Perhaps more importantly, Expedition Fiord is well mapped compared to most of the CAA, where few previous small-scale topographic products exist to compare with recent data. An ability to directly use 1950/60s aerial photography in current SfM-MVS processing would provide more accurate historical topographic reconstructions, and allow geodetic mass-balance estimates to be made for previously unstudied glaciers (Mölg and Bolch, 2017). This method has been used to quantify long-term change on individual glaciers in Central Asia (Denzinger and others, 2021) and the Swiss Alps (Mölg and Bolch, 2017), as well as on a regional scale in the Argentinian Andes (Falaschi and others, 2022) and on Svalbard (Geyman and others, 2022), but its potential remains largely underutilised in the CAA.

Here, we assess long-term changes of Bowman Glacier, a small mountain glacier located in Tanquary Fiord, northern Ellesmere Island (Fig. 1). We use SfM-MVS processing applied to recent and historical aerial photography to reconstruct glacier surface topography in 1959 and 2018, producing detailed point clouds, DEMs and orthomosaics for the glacier. These provide information concerning long-term ice surface elevation changes and allow determination of the geodetic mass balance of the glacier over a 60 year period. This is combined with optical satellite data to reconstruct the evolution in extent of the glacier over the past six decades, and ground-penetrating radar (GPR) measurements of ice thickness to estimate the remaining ice volume. This information is compared with Regional Atmospheric Climate Model (RACMO) surface mass-balance data from the last six decades, and changes are projected into the future to estimate the life expectancy of the remaining ice mass.

Fig. 1. Study area near Tanquary Fiord, Ellesmere Island, with location of Bowman Glacier (dashed white box) and its three main drainage valleys (white lines), including May Creek which provides the main water supply for Tanquary Camp. Base image: Sentinel 2A, 24 July 2020.

Study site

Bowman Glacier (81.35$^\circ$N, 76.45$^\circ$W) currently consists of several small disconnected ice masses, located ~10 km southeast of the Parks Canada basecamp at the head of Tanquary Fiord, Ellesmere Island, Nunavut, Canada (Fig. 1). Established in 1963 by the Defence Research Board as a scientific field station, the camp now serves as the headquarters of Quttinirpaaq National Park. Lying on a plateau that reaches ~1200 m above sea level (a.s.l.), Bowman Glacier is located on the western edge of the Hazen Plateau, a sparsely glaciated area between the Northern Ellesmere Icefields to the north, and Agassiz Ice Cap to the south. When Tanquary Camp was established, Bowman Glacier drained into three main valleys, with the eastern and central drainage basins feeding into Macdonald River to the north, and the western one providing the main water supply for the camp via May Creek (Figs 1, 2). The initial motivation for this study was to estimate the remaining life expectancy of Bowman Glacier and to determine its viability as the primary water source for Tanquary Camp.

Fig. 2. (a) Oblique air photo (T409R-178) from 16 July 1950 (Hattersley-Smith, 1969), with location of Bowman Glacier, May Creek and Tanquary Camp. Lake Hazen, ~80 km to the northeast, is visible in the top left corner. (b) Oblique air photo taken during the helicopter survey on 3 August 2018. Dashed white box in (a) and (b) indicates the same region as in Fig. 1.

Previous observations on the glacier include only a ‘limited assessment of snow cover and ablation’ from 1966 (Hattersley-Smith and Oliver, 1967). Although Hattersley-Smith (1969) briefly mentions Bowman Glacier as part of a regional summary of glacial features, there is no previous published work on the ice mass. The earliest known images of the glacier are from a trimetrogon survey flown by the Royal Canadian Air Force in July 1950 (Cogley and Adams, 2000), which include both nadir and oblique photos (Hattersley-Smith, 1969) (Fig. 2a). The primary photographic record of Bowman Glacier was acquired in summer 1959, when the Royal Canadian Air Force completed comprehensive nadir photography flights across the CAA at a scale of 1:60 000, which were used to derive the standard 1:250 000 and 1:50 000 scale maps of the region.

Materials and methods

Ice extent

A combination of aerial photographs and satellite imagery was used to reconstruct the ice extent of Bowman Glacier between 1959 and 2020 (Table 1; Fig. 3). The 1959 extent was manually outlined from an 8 m panchromatic orthomosaic (Fig. 3a) generated from air photos using SfM-MVS processing as detailed further in the following section. A 10 m resolution SPOT-1 (Satellite pour l'observation de la Terre) panchromatic image was used to manually outline the 1987 extent (Fig. 3b; Table 1). All manual digitisation was performed in QGIS (version 3.12.2).

Fig. 3. Evolution of Bowman Glacier over the period 1959–2020 from aerial photography (a), and optical satellite imagery (b–g). Panchromatic images: (a) orthomosaic from August 1959 air photos; (b) SPOT-1 20.08.1987; false-colour composites: (c) Landsat 7 24.06.2000; (d) Landsat 7 25.07.2005; (e) Landsat 7 30.06.2010; (f) Landsat 8 25.07.2015; (g) Landsat 8 17.08.2020.

Table 1. Optical satellite imagery used for mapping ice extent

SPOT-1 HRV: Satellite pour l'observation de la Terre High resolution visible; ETM+: Enhanced Thematic Mapper Plus; OLI: operational land imager; SWIR: shortwave infrared.

Changes in ice extent since 2000 were mapped at 5 year intervals from Landsat 7 (2000–10) and Landsat 8 (2015–20) multispectral imagery (Table 1) using raw digital number data and the normalised difference snow index (NDSI). The NDSI takes advantage of the contrast in spectral reflectance between snow-covered terrain and surrounding rock areas, and is calculated from the green and shortwave infrared (SWIR) bands from: (green−SWIR)/(green+SWIR) (e.g. Paul and others, 2015). Orthorectified Landsat data (Level 1T) were accessed through the Google Earth Engine service and processed using a custom script and the online JavaScript code editor. For each year, all suitable cloud-free images acquired in the summer months (June–August) were reprojected to a common coordinate reference system (WGS 84, UTM 18N) and used to compute band ratio images, which were then combined to generate a single median NDSI raster with 30 m pixel size. An assigned threshold value separating glacier ice from surrounding terrain was used to derive binary images indicating the presence or absence of ice.

For 2015 and 2020, there were >25 suitable Landsat 8 scenes available for each summer season, compared to <10 Landsat 7 scenes for earlier years. Using a multitemporal image stack helps fill data gaps from the banding caused by the scanline corrector failure on the ETM+ instrument starting in 2003. Another advantage of using multiple images is a lower sensitivity to variations in illumination such as shadows due to changing sun angles, and differences in snow cover due to persistent snowpatches in early summer and snowfall events later in the season. This is especially effective for 2015 and 2020 with the increased availability of Landsat 8, as the median NDSI remains more stable between years and the threshold value separating the glacier surface from ice-free terrain was kept constant at 0.10, with any pixel with a value above the threshold being classified as ice. With fewer images available for 2000–10, median NDSI values were more variable so the threshold was adjusted between 0.28 and 0.32. Additional factors accounting for the difference in assigned thresholds for OLI vs ETM+ imagery include the difference in the wavelengths of the green and SWIR bands (Table 1), and different gain values used for the two sensors.

In order to remove isolated pixels and smooth object boundaries, a morphological opening operation, consisting of an erosion followed by a dilation, was performed on the binary images. With both operations applied within the same three-pixel kernel, the erosion shrinks clusters of contiguous pixels and removes small objects, while the dilation restores larger objects, mostly preserving their shape and size. It is effective at reducing noise due to patches of seasonal snow being incorrectly classified as ice. The smoothed binary images were then vectorised to produce glacier outlines, which were exported from Google Earth Engine for further processing in QGIS.

On years when significant snow cover persisted throughout the summer (i.e. 2000 and 2005; Figs 3c, d), the opening operation failed to filter some of the larger snowpatches. In those cases, we used the 1959 glacier outlines, corresponding to the maximum ice extent for the study period, and excluded all individual patches of pixels falling entirely outside this boundary. Like other automated classification methods, the NDSI struggles with correctly identifying glacier ice in strong shadows, and some outlines (particularly near the terminus) had to be manually adjusted based on visual inspection of false-colour composite Landsat images. Additional high-resolution Digital Globe imagery from the summers of 2008 (Quickbird 2, 8 June 2008, 0.6 m resolution) and 2011 (Worldview 2, 22 July 2011, 0.5 m resolution) was used as reference to help differentiate perennial snow from glacier ice. We calculate uncertainties in glacier area using a ±1/2 pixel buffer, similar to Granshaw and Fountain (2006). Relative errors increase with decreasing glacier size, and range between 4% of glacier extent for the higher resolution 1987 SPOT-1 image and 16–24% for the coarser 30 m Landsat data for 2000–20.

Ice surface elevation

The full SfM-MVS workflow summarised below was performed in Agisoft Metashape Pro (version 1.6). Additional point cloud filtering and distance computations between the 1959 and 2018 point clouds were undertaken in CloudCompare (version 2.11.1). The final point clouds were used to derive georeferenced DEMs and orthomosaics of the study area.

Ice surface reconstruction: 1959

The 1959 surface topography was determined from a series of panchromatic 1:60 000 scale vertical aerial photographs archived at the National Air Photo Library in Ottawa, Canada. Images were acquired on 29 July 1959 (roll A16691, photos 49–52) and 26 August 1959 (roll A16788, photos 62–65 and 77–80), from an altitude of 9000 m a.s.l. using a Wild RC5A aerial film camera and Aviogon lens (15 AG) with a nominal focal length of 153 mm. The 9 in × 9 in (22.9 cm × 22.9 cm) images were scanned 20 years ago, using a desktop flatbed scanner at 300 dpi, resulting in a pixel pitch of 85 μm and an equivalent ground sampling distance of ~5 m.

The scanning process is a limiting factor in this study, with the main drawbacks being relatively low spatial resolution (300 dpi) and geometric accuracy of the scanned material. The recommended minimum scanning resolution for aerial photography is 600–1200 dpi, and high accuracy photogrammetric applications require photogrammetric-quality scanners with sub-pixel geometric accuracy of 2–5 μm. Many desktop scanners, such as the one used in this study, have a typical geometric accuracy of 50–100 μm (Baltsavias and Waegli, 1996) which, at a 1:60 000 scale, gives planimetric errors of 3–6 m, roughly equivalent to 1 pixel on the ground (5 m). This is sufficient for the requirements of this study, but could be improved upon for future work.

A total of 12 overlapping air photos from three different flightlines (two of which were flown on the same day) were imported into Metashape, where fiducial marks (four marks, one in each corner) were manually placed on each photo and used to automatically align and scale photos to the same dimensions. A single mask was applied to the full dataset to exclude the outer frame from the scanned photos.

After the image alignment step, the sparse cloud was georeferenced to the 2018 DEM based on 13 control points in proximity to the 1959 ice extent (Fig. 4a). Distinct stable bedrock features were manually selected in the 2018 dense point cloud, and subsequently identified in the 1959 air photos. Features were selected to ensure positional accuracy of 10 m or less (the equivalent of ±2 pixels), and each was positioned on 4–6 photos. The georeferenced cloud was optimised and points with reprojection errors >0.25 pixels were deleted (~10% of original point count). The final cloud contained ~49 000 points with RMS reprojection errors of 0.24 pixels. Combining all 13 control points, RMS registration errors were 6.53 m in total, or 5.75 m in horizontal and 3.09 m in vertical. The dense point cloud, generated with the ultra-high quality setting and aggressive filtering in Metashape, contained nearly 16 million points.

Fig. 4. (a) Study area with flightpath of the 3 August 2018 air photo survey and locations of the two ground targets (numbered 1 and 2) used as checkpoints for the direct georeferencing method, and the 13 control points (white dots) used to align the 1959 and 2018 reconstructions. Base image: Sentinel 2B, 3 August 2018. (b) Image captured during the air photo survey (3 August 2018), with insets of surveyed targets (1 m × 1 m checkered flags).

Ice surface reconstruction: 2018

The aerial photo survey took place on 3 August 2018 onboard an Airbus ASTAR 350 B2 helicopter. We surveyed an area of ~70 km² over an elevation range of nearly 1000 m, including the Bowman Glacier summit plateau and the three main valleys in its drainage basin. The aircraft maintained a constant altitude of ~1400 m a.s.l., or between 50 and 1350 m above ground level, covering a linear distance of ~85 km in just under an hour of flying time. Over 1000 oblique images were captured with a Nikon D850, a full-frame (35 mm) digital single-lens reflex camera, and a 24 mm prime lens (NIKKOR AF-S 24 mm f/1.8G ED).

Georeferencing was achieved using direct control measurements acquired with a GPS, including a survey-grade dual frequency (L1/L2 bands) Trimble R7 receiver and Zephyr 1 Geodetic antenna. Prior to the scheduled survey, the GPS was used to record the precise position of two checkered flags (1 m × 1 m) with a minimum occupation time of 30 min (Fig. 4b). The two targets were used as checkpoints for validating the direct georeferencing method. The GPS was initially set to record satellite observations throughout the flight at a 10 Hz logging rate, but due to a system malfunction the receiver reset itself to factory defaults and only recorded observations at 15 s intervals. Synchronisation between the camera shutter and positioning system was achieved via the Trimble Event Input Marker, which connected the GPS receiver to the camera hot shoe accessory mount. Camera positions were subsequently calculated using interpolation based on the moment of image capture recorded by the GPS receiver.

Optimised image files and corresponding masks were imported into Metashape and used to generate a sparse point cloud, which was subsequently georeferenced using the camera position information assigned to each image. About 10% of tiepoints with the highest reprojection error values were removed and, following optimisation, RMS point reprojection errors were 0.41 pixels. Local RMS registration errors, based on the two checkpoints measured in the field, were 0.35 m horizontal and 0.68 m vertical (0.77 m total). The dense point cloud was generated in Metashape using the ‘high’-quality setting and the mild filtering option, in order to reduce processing time and remove some of the biggest outliers. Compared to the 6 million points in the sparse cloud, the final reconstruction contained 1.56 billion points.

Surface elevation change: 1959–2018

Point cloud data generated from image-based reconstruction techniques are inherently noisy, with outliers and geometrical errors arising from imperfections in the source data (e.g. noisy or out of focus images) or inadequate scene coverage. Filtering this surface noise and removing outliers helps to reduce artefacts and enable geometrically consistent surface reconstructions. In this study, filtering was performed using the Multiscale Model to Model Cloud Comparison (M3C2; Lague and others, 2013) plugin implemented in CloudCompare. Developed for direct cloud comparison and change detection, the M3C2 algorithm can also be used to filter point cloud data by computing local statistics. The algorithm calculates the median height values of all points within a cylindrical neighbourhood of a given diameter and depth, oriented along a projected vector perpendicular to the surface (i.e. surface normal).

Filtering parameters in M3C2 were selected to take into account surface roughness and point density of the reconstructions, and to ensure a minimum of five points within each neighbourhood. For the 1959 point cloud, with a mean point spacing of ~8.3 m and >99% of points being 10 m or less from their nearest neighbour, surface normals were computed over a scale of 100 m, the projection diameter was set at 20 m and maximum depth at 50 m. Given the mean point spacing, the final DEM and panchromatic orthomosaic were gridded at 8 m resolution, with each pixel representing the mean of all elevation or intensity values of all points within a given gridcell. For the 2018 point cloud, surface normals were calculated at varying scales, between 2 and 10 m. Given a mean point spacing of ~0.45 m in 2018, the projection scale was set at 1 m and the depth at 5 m. The smoothed point cloud was then gridded into DEM and RGB orthomosaic rasters at 0.5 m horizontal resolution (Fig. 5).

Fig. 5. SfM-MVS workflow outputs from the 3 August 2018 air photo survey. (a) Elevation model of the study area including Bowman Glacier (white outline) and its drainage system with stream network (blue lines) and drainage basin boundaries (white lines). Contour interval of 20 m. (b) RGB orthomosaic of the study area.

Following the filtering step, both the 1959 and 2018 point clouds were cropped to the same extent shown in Fig. 6, covering an area of 17 km² including Bowman Glacier and the summit plateau. Differences in ice surface elevation between 1959 and 2018 were then calculated by comparing the two point clouds in 3D using the M3C2 algorithm (Lague and others, 2013). As opposed to DEM differencing, performing change detection directly on point clouds allows the gridding step to be skipped, and therefore avoids the introduction of additional uncertainties from interpolation of elevation values. Surface change was estimated from median point elevations calculated along the same surface normals computed on the 1959 cloud in the filtering step, while also keeping the same projection scale and maximum search depth (20 and 50 m).

Fig. 6. (a) Panchromatic orthomosaic of Bowman Glacier generated from the 1959 air photos, showing the 1959 and 2018 glacier extents (black outlines). (b) Ice surface elevation change between 1959 and 2018 from dense point cloud comparison. Isolines of ice surface elevation change at 5 m. Contour interval over stable terrain of 10 m.

Uncertainties were determined taking into account a 3 m registration error, corresponding to the vertical component of the misalignment errors calculated over the 13 control points used to align the two clouds. In combination with this 3 m error applied uniformly across the entire cloud, M3C2 uses local surface roughness characteristics computed at the projection scale to calculate a spatially variable level of detection at 95% confidence for each distance measurement, which can be used to determine the statistical significance of any detected change (Lague and others, 2013). The corresponding data were gridded into 8 m resolution rasters which were then used to calculate ice thickness and volume loss estimates.

Ice thickness from GPR measurements

The current volume of Bowman Glacier was estimated from spot ice-thickness measurements made with a 10 MHz GPR system on 3 August 2018. The system consists of a monopulse transmitter based on the design of Narod and Clarke (1994), and a TekTronix THS720A digital oscilloscope as the receiver, with images of the radar returns recorded with a digital camera. A handheld GPS unit was used for horizontal positioning, with an estimated accuracy of <10 m; the elevation of the measurement points was extracted from the 2018 SfM DEM. A total of 11 points, spaced ~100 m apart, were measured along two transects, one running east-west across the top of the main ice mass, the other north-south along the centreline. A radio-wave velocity of 0.171 m ns⁻¹ was used to convert GPR returns to true ice thickness, based on the mean of 12 common mid-point surveys undertaken by Copland and Sharp (2000) at John Evans Glacier, 200 km south of Bowman Glacier.

Surface mass balance

Surface mass balance (SMB) was derived from the regional climate model (RACMO) version 2.3 following the method detailed by Noël and others (2018). To more accurately resolve complex topography and small isolated ice masses such as Bowman Glacier, daily SMB components (precipitation, melt, runoff, refreezing, sublimation and drifting snow erosion) were statistically downscaled from the original 11 km RACMO2.3 model output to 1 km resolution. Elevation and albedo biases affecting surface melt and runoff components were corrected using the Canadian Digital Elevation Model (CDEM), glacier outlines from the Randolph Glacier Inventory 5.0 (Pfeffer and others, 2014), and bare ice albedo information derived from MODerate-resolution Imaging Spectroradiometer (MODIS) data. Model results were evaluated and calibrated using field measurements from 1959 to 2010 from 198 sites across the northern CAA. Daily SMB estimates were obtained by summing the individual components and combined into monthly values spanning the 1958–2019 period. In this dataset, Bowman Glacier is represented by a single 1 km² pixel approximately centred over the main basin, and any changes in ice extent and glacier hypsometry due to mass loss over the study period are therefore not accounted for.

Results

Ice extent and thickness change

Ice cover mapping from aerial photography and satellite imagery shows that Bowman Glacier shrank over the past six decades (Figs 3, 7). In 1959, the glacier consisted of seven independent ice masses, of which three (including the main one) drained into May Creek. By 2000, the two westernmost ice masses had virtually disappeared, and only a small portion (~0.058 km²) of the main glacier remained in the May Creek basin. By 2005 this remnant in May Creek had also disappeared, leaving glacier ice in only the central and eastern basins, which drain to Macdonald River and not to Tanquary Camp. Over the past 61 years, the total ice-covered area has decreased from 2.75 km² in 1959 to 0.61 km² in 2020, corresponding to a 78% area loss. In 41 years between 1959 and 2000, total ice cover loss was 0.89 km², while in the period 2000–20 total area loss was 1.26 km², meaning that average annual retreat rates nearly tripled between the two periods (0.022 vs 0.063 km² a⁻¹).

Fig. 7. Evolution of ice extent of Bowman Glacier for 1959–2020. The white lines show the boundaries of the three drainage basins, with the western one providing the main water supply to Tanquary Camp. In 2000 a small remnant ice patch still drained into May Creek. By 2005 all the ice had retreated out of the western basin. Contour interval of 10 m.

Results from differencing the 1959 and 2018 reconstructed ice surface show average elevation changes over the glacier of −14.7 m with a std dev. (σ) of 9.2 m, and maximum losses reaching >40 m near the main terminus (Fig. 6b). Losses also reached >30 m in the lower regions of the easternmost ice masses. Random error in the elevation change measurements is calculated from apparent change computed over an 11.35 km² area of stable ground adjacent to Bowman Glacier, including all pixels within a 1 km distance of the glacier limits. The average elevation difference over this whole area is −0.03 m, and the RMS error 2.8 m. Figure 8 shows the hypsometry for 1959 and 2018 in 25 m elevation bands, as well as ice thickness and volume change for the same bands. In terms of ice cover, the largest losses occurred in the middle elevation bands of 1075–1125 m a.s.l., where $\lt$20% of the 1959 ice-covered area remained in 2018. As most of the 1959 ice cover was concentrated above 1100 m, total volume losses are highest in those elevation bands. In comparison, surface lowering is highest between 1000 and 1050 m a.s.l., with average elevation losses reaching >20 m (Fig. 8).

Fig. 8. (a) Bowman Glacier hypsometry in 25 m elevation bands for 1959 (light grey bars) and 2018 (dark grey bars), and average ice surface elevation change (blue line) for the same bands. (b) Average ice-volume loss for 1959–2018 per elevation band. Values on the vertical axis represent the midpoint elevation of each band.

GPR measurements from 2018 show an average ice thickness of ~25 m in the upper part of the basin (previously the accumulation area of the glacier), and ~42 m in the lower part, with an overall average of 34.1 m. Comparison of the 1959 and 2018 ice surface elevations along the GPR transects suggests that the mean ice thickness in 1959 would have been 44 m in the upper part, and 65 m lower down, for an average of 55.5 m. This corresponds to an average ice surface elevation loss of 21.5 m, meaning that the mean ice depth remaining in 2018 was $\sim$60% of the 1959 thickness (Fig. 9).

Fig. 9. (a) Location of 11 GPR spot measurements performed along two transects on 3 August 2018. (b, c) Ice surface elevation change between 1959 (light grey) and 2018 (dark grey) along two GPR transects: north-south along the glacier centreline (b), and east-west across the top of the main ice mass (c). The black dashed line represents the approximate bed elevation along each transect interpolated between ice-thickness measurements (black dots).

Ice-volume and mass loss

The 1959–2018 total ice-volume loss over the full extent of the glacier, determined as the sum of the product of height differences for all pixels within the 1959 glacier extent and their corresponding area, is 0.041 ± 0.0083 km³. The uncertainty term (±0.0083 km³) is based on 95% confidence interval calculated for each cloud-to-cloud distance measurement within the 1959 extent. Given the absence of snow cover on the ice surface on the 1959 air photos and in more recent satellite imagery from 2010 to 2020, it is reasonable to assume that the full extent of Bowman Glacier lies in the ablation area, and consequently has no significant accumulation zone. We therefore convert volume to mass change assuming a constant density profile and using a conversion factor approaching ice density, of 900 kg m⁻³ (Huss, 2013), and calculate a total loss of 0.036 ± 0.008 Gt. To compensate for glacier retreat over the 59 year period, we take the average ice area calculated for 1959 and 2018 (1.78 km²) and determine a glacier-wide average thickness change of −22.7 ± 4.7 m. This corresponds to a total of −20.5 ± 4.2 m w.e., and an average of −347.0 ± 71.4 mm w.e. a⁻¹ over the 1959–2018 period.

Surface mass balance

Modelled monthly surface mass-balance data for 1958–2019 show a cumulative loss of 24.56 m w.e., with a noted increase in the rate of mass loss in the last two decades (Fig. 10). Average losses more than doubled between 1958–99 and 2000–19, increasing from 279 mm w.e. a⁻¹ (σ = 329 mm w.e. a⁻¹) in the first period, to 642 mm w.e. a⁻¹ (σ = 401 mm w.e. a⁻¹) in the second (Table 2). While most years have been negative, the frequency of positive net balance years decreased from an average of 2.5 per decade prior to 2000, to none since then, and four out of the five most negative years occurred in the last decade. With a few exceptions, there is overall little variability in winter balance over the study period (σ = 40 mm w.e. a⁻¹), with a fairly constant accumulation pattern and mean winter balance values decreasing slightly after 2000 (~20% decrease, from 142 to 117 mm w.e. a⁻¹). In contrast, there is a marked variability in summer balance (σ = 365 mm w.e. a⁻¹) associated with enhanced melt rates in warmer summers. Average summer losses increased from 422 mm w.e. a⁻¹ in 1958–99 to 759 mm w.e. a⁻¹ in 2000–19, corresponding to an ~80% increase since the start of the 21st century (Fig. 10).

Fig. 10. Winter (white), summer (light grey) and net annual (dark grey) SMB for 1958–2019, with 10 year running mean (dotted line) and decadal averages (thick horizontal lines). Cumulative total SMB (solid line) is shown on the secondary axis. SMB data derived from the regional climate model RACMO2.3 downscaled to 1 km resolution, and values determined over a 1 km² gridcell centred over Bowman Glacier.

Table 2. Net, winter, and summer SMB for Bowman Glacier over the 1958–2019 period, derived from the regional climate model RACMO2.3

Decadal and multi-decadal averages ($\overline {x}$) and std dev. (σ) in units of mm w.e. a⁻¹.

Discussion

Between 1959 and 2020 Bowman Glacier shrank by 2.15 km², losing 78% of its original extent, leaving it with an area of only 0.61 km². The majority of this loss, ~1.26 km², occurred in the last two decades, at a rate almost three times higher than before 2000. Overall, there are few studies on small ice masses in the northern CAA, although there are mass-balance measurements from two pairs of small ice caps on the Hazen Plateau, Murray and Simmons and St. Patrick Bay ice caps, 100–200 km east of Bowman Glacier (Bradley and Serreze, 1987; Braun and others, 2004; Serreze and others, 2017). In 2016, the four ice caps had lost between 61 and 95% of their 1959 extents (2.94–7.45 km²). With the smaller of the two then only covering 0.15 km², Serreze and others (2017) predicted that both St. Patrick Bay ice caps would be gone within 5 years; satellite imagery from 2020 indeed shows both glaciers have now completely disappeared. Using remote-sensing data, Curley and others (2021) reconstructed the evolution of eight land-terminating glaciers in Alexandra Fiord, 280 km south of Bowman Glacier on the east coast of Ellesmere Island, and found average annual area loss rates to have doubled between the 1959–2001 and 2001–19 periods. The smallest of the eight glaciers shrank by ~50% since 1959 (from 1.81 to 0.88 km²), compared to a total area loss of ~12% over the same period when including all eight glaciers (of which seven were >5 km² in 2019).

On a regional scale, small ice masses have been experiencing the highest relative area losses, at least for the last six decades (Sharp and others, 2014; White and Copland, 2018). In the only recent assessment of regional changes in glacier ice extent over Northern Ellesmere Island, White and Copland (2018) found that, between 1999 and 2015, the rate of area change for small land-terminating glaciers <1 km² in size averaged $-$30.4% decade⁻¹, compared to −0.4% decade⁻¹ for larger glaciers >1000 km². In comparison, between 1959/60 and 2000/01, ice masses <1 km² decreased in area by an average of 29.2% (7.3% decade⁻¹) and those in the size class above (1–5 km²) by 25.6% (6.4% decade⁻¹) (Sharp and others, 2014). Averaged over the same periods, the rate of area change on Bowman Glacier increased from −7.9% decade⁻¹ between 1959 and 2000 to −36.3% decade⁻¹ between 2000 and 2015, which is comparable to the regional average for ice masses in similar size classes (Sharp and others, 2014; White and Copland, 2018). Besides glacier size, White and Copland (2018) also identified elevation and glacier flow length as the main controls on area changes, with smaller and shorter glaciers at lower altitudes experiencing highest relative area losses. This is similar to findings from Alexandra Fiord where, among the eight studied glaciers, the two with the smallest elevation range, and with most of their area concentrated at lower elevations, lost the greatest proportion of their original extent, shrinking by $\sim$51% and $\sim$34% between 1959 and 2019 (Curley and others, 2021).

Thinning and mass loss

Our results for Bowman Glacier show a glacier-wide ice surface lowering of 22.7 ± 4.7 m between 1959 and 2018, giving average thinning rates of 0.39 ± 0.08 m a⁻¹. This is comparable to results from the only available regional assessment of ice-thickness change in the northern CAA, where Hugonnet (2021) reported average thinning rates on land-terminating glaciers of 0.39 ± 0.03 m a⁻¹ for the 2000–19 period, including a 38% increase in the rate of change between 2000–04 and 2015–19. Although comparable, the surface lowering rates determined for Bowman Glacier are averaged over a considerably longer period (1959–2018), which includes a few years of positive mass balance in the 1960s to 1990s (Fig. 10). Considering only the last two decades, modelled SMB on Bowman Glacier shows average surface lowering of 0.71 m a⁻¹, and a 90% increase in rate of thinning between 2000–04 and 2015–19, meaning that both overall annual losses and interannual thinning rates are greater than the regional average for land-terminating glaciers determined in Hugonnet (2021).

Differencing the 1959 and 2018 reconstructed surfaces for Bowman Glacier yields a total volume loss of 0.036 ± 0.008 Gt, and a glacier-wide average elevation change of −20.5 ± 4.2 m w.e., at a mean specific rate of −347.0 ± 71.4 mm w.e. a⁻¹. In comparison, modelling results from RACMO2.3 averaged over the same period (September 1959 to July 2018 inclusively) show a net SMB of −22.0 m w.e. and a mean specific rate of −372.8 mm w.e. a⁻¹, equating to a 7.4% difference from our estimates, well within the error margin. Total volume losses for the 1 km RACMO2.3 gridcell sum to 0.022 Gt. When scaled to the mean area of Bowman Glacier of 1.78 km², this amounts to 0.039 Gt, which is in close agreement with the 0.036 ± 0.008 Gt total mass loss determined from differencing the 1959 and 2018 point clouds.

The reconstructed hypsometry at Bowman Glacier for 1959 and 2018 shows net area losses and ice surface lowering at all elevations (Fig. 8). Geodetic mass-balance measurements for 1960–2014 at White Glacier, a 14 km long valley glacier on western Axel Heiberg Island, 340 km south-west of Bowman Glacier, show overall thinning at elevations below 1400 m a.s.l., with maximum losses of ~40 m near the terminus (Thomson and others, 2017). Between 900 and 1300 m a.s.l., roughly the elevation range of Bowman Glacier, ice-thickness losses at White Glacier ranged between 10 m and 23 m over the period 1960–2014, which is comparable to the thinning values computed for the equivalent elevation bands at Bowman Glacier of 13–21 m (Fig. 8).

At White Glacier, the mean geodetic mass balance for 1960–2014 averaged over a glaciated area of 39.81 km² is −9.61 m w.e., or −178 ± 16 mm w.e. a⁻¹ (Thomson and others, 2017), about half the losses determined at Bowman Glacier in this study. There are notable differences between the two glaciers, including their length, elevation span and the fact that, unlike Bowman Glacier, White Glacier still has an accumulation area which offsets some of the ice loss observed at lower elevations. Bowman Glacier now measures <1 km along its main axis, covers an elevation range of only 350 m (between 900 and 1250 m a.s.l.), and is entirely in the ablation area. In contrast, White Glacier is 14 km long and spans an elevation range of 1700 m, with the terminus descending to nearly sea level and a 5 km wide accumulation area reaching 1780 m a.s.l. (Thomson and others, 2017).

Connections between mass loss, hypsometry and regional equilibrium line altitude

Aerial photos from late August 1959 (Figs 3a, 6a) show bare ice over the entire extent of Bowman Glacier, with the snowline only visible on the much larger ice cap 2 km to the south, at an altitude of ~1250 m. Regional reconstructions of the lowest equilibrium line altitude (ELA) over the northern CAA in 1960 show it reaching 1200 m over the northern Ellesmere Icefields to the north of Tanquary Fiord, and 1000–1100 m over the Hazen Plateau to the east (Miller and others, 1975; Wolken and others, 2008), which would place the summit of Bowman Glacier either at the limit, or just below, the local ELA at that time.

With the exception of a period of cooler conditions in the 1970s, a number of studies point to a general increase in ELA since the 1960s across the northern CAA. For example, on the eight glaciers in Alexandra Fiord, the mean end of summer snowline elevation (assumed to represent the ELA) increased by 77 m between 1974–82 and 2014–19 (Curley and others, 2021). Field measurements at White Glacier and on Meighen, Agassiz, and Devon ice caps, show the average decadal ELA to have increased by ~330 m between the 1980s and 2010s (Thomson and others, 2017; Burgess and Danielson, 2022). White and Copland (2018) demonstrated the impact of rising ELA on the ice masses of northern Ellesmere Island and showed that a 100 m rise above the 1960 ELA would cause the regional ablation area to increase from 35.2 to 47.1% of total glacier cover. For a 300 m rise, the ratio would increase to 72.3 and >50% of all glaciers in the region would be entirely below the ELA (White and Copland, 2018). Glaciers at relatively low elevations, and those in areas with relatively flat topography with the bulk of ice within a narrow elevation band (e.g. small glaciers on plateaus, such as Bowman Glacier), are more sensitive to changes in regional climate, particularly when the ELA rises above their summit elevation. For example, of the four small ice caps on the Hazen Plateau, the two that recently disappeared were located 200–300 m lower in elevation than the nearby Murray and Simmons ice caps, which still survive (Bradley and Serreze, 1987; Braun and others, 2004; Serreze and others, 2017). However, similar to Bowman Glacier, Murray and Simmons ice caps were already at the limit of the local ELA estimated at 1000–1100 m a.s.l. in 1960 (Miller and others, 1975; Wolken and others, 2008), and so their loss is inevitable under the current climate. The rapid shrinkage of small glaciers at lower altitudes in the study region (Sharp and others, 2014; Serreze and others, 2017; White and Copland, 2018) highlights the importance of glacier hypsometry in connection to regional ELA, and its role in determining mass-balance sensitivity to changes in climatic conditions.

Changes in regional climate

The regional climate in the northern CAA is characterised by relatively low annual precipitation (<400 mm a⁻¹), with some variability depending on the combined effect of topography and distance to a moisture source (Koerner, 1979). Precipitation is highest along coastlines facing Baffin Bay and the Arctic Ocean, and decreases with elevation and distance inland. These precipitation gradients in part explain the difference in local ELA between Alexandra Fiord (~670 m a.s.l.), White Glacier (~1075 m a.s.l.), and Bowman Glacier (~1200 m a.s.l.). The 1950–2020 instrumental record at the two nearest weather stations, Eureka and Alert, 200 km southwest and 250 km northeast of Bowman Glacier, respectively, shows average annual air temperatures of –18.9$^\circ$ and –17.4$^\circ$C, and mean total annual precipitation of 72 and 156 mm (Environment and Climate Change Canada, 2022). Both stations recorded an average warming of >1.5$^\circ$C between 1950–99 and 2000–20 (Fig. 11a), with winter and autumn temperatures showing the largest increase (>2$^\circ$C) between the two periods (Fig. 12b).

Fig. 11. (a) Mean annual air temperature and (b) total annual precipitation for Eureka (black) and Alert (grey) weather stations. Annual values (dotted lines) and 10 year running mean (solid lines) for 1950–2020. Data source: Environment and Climate Change Canada (2022).

Fig. 12. (a) Mean summer and (b) winter air temperature for Eureka (black) and Alert (grey) weather stations. Annual values (dotted lines) and 10 year running mean (solid lines) for 1950–2020. Data source: Environment and Climate Change Canada (2022).

The annual mass balance of glaciers and ice caps in the CAA mainly depends on summer melt rates rather than changes in winter accumulation (Koerner, 2005), and the accelerated mass loss in the last two decades has been shown to coincide with an increase in summer (June, July, August) air temperatures (Sharp and others, 2011; Mortimer and others, 2016). At Bowman Glacier, there is a strong positive correlation (R ² = 0.99, p − value = 2.22 × 10⁻¹⁶) between yearly summer losses and net balance, indicating that the increasingly negative mass budget is in direct response to enhanced summer melt due to warmer summer temperatures. At Eureka, total annual precipitation has increased by $\sim$30% since 1950 ($\sim$4% decade⁻¹), but there is no significant trend in the precipitation record at Alert (Fig. 11b). There is little interannual variability in winter mass balance over Bowman Glacier (Fig. 10), suggesting that any potential increase in total precipitation is not enough to offset losses from warming temperatures.

A similar pattern of enhanced mass loss has also been observed on peripheral glaciers in North Greenland, where air surface temperatures have been rising faster than further south, and thinning is now occurring at all elevations. Particularly warm summer conditions in recent years, due to enhanced atmospheric advection from the mid-latitudes, have contributed to a fourfold increase in glacier mass loss over the 2003–21 period, and thinning rates in North Greenland are now higher than for all other peripheral regions combined (Khan, 2022). Melt rates are expected to increase across the Arctic in the following decades due to continued warming and lengthening of the melt season (IPCC, 2022). Albedo feedback processes will enhance these losses, particularly on small glaciers, due to the exposure of dirt layers and darkening of the ice surface and to increased heat advection from exposure of surrounding bedrock as the ice cover shrinks (Sharp and others, 2011; Serreze and others, 2017). Models predict rising surface air temperatures accompanied by increased mean precipitation, but with a higher fraction falling as rain (IPCC, 2022), this also contributes to ablation and lowers albedo.

Future projections

Given that Bowman Glacier has lost nearly 80% of its area over the past 60 years, and that the entire glacier now sits within the ablation zone, the ice mass is expected to completely disappear over the coming decades. The predicted timing of its demise varies depending on the parameter used as a basis for the calculations, namely:

1. (1) Volume change: taking the 2018 glacier area (0.82 km²) and the average measured ice thickness along the two GPR transects (34.1 m), gives a total remaining ice volume of 0.028 km³, or 0.025 Gt. Based on our estimates, the glacier lost 0.041 ± 0.008 km³, or ~60% of its volume since 1959. Projecting the average rate of volume change determined for 1959–2018 (−0.0069 km³ decade⁻¹) into the future suggests that the glacier will completely disappear in just under four decades, by ~2060.

2. (2) SMB: from the modelled SMB data, the average thinning rate for the last decade was 709.2 mm w.e. a⁻¹. Taking the average ice thickness measured in 2018 (34.1 m of ice, 30.69 m w.e.), and assuming that this rate continues into the future, suggests that the disappearance might occur by ~2060.

3. (3) Area change: in 2020, the last year of our study period, Bowman Glacier covered only 0.61 km². Projecting the average rate of change in ice extent from the last decade (−0.057 km² a⁻¹) beyond 2020, suggests Bowman Glacier will have completely disappeared by ~2030 (Fig. 13).

Fig. 13. Change in ice extent (black), volume (dark grey) and average thickness (light grey) of Bowman Glacier for 1959–2020 (solid lines) and projected rate of change (dotted lines), represented as remaining percentage in reference to 1959. Volume projections are based on the 1959–2018 average, while area and thickness change are based on the 2010–20 trend.

Bowman Glacier is committed to continued mass loss, regardless of climate scenarios. The above projections likely bound the lower and upper limits for the complete disappearance of the glacier, sometime between 2030 and 2060, provided that regional warming trends stabilise. However, given that the most recent Intergovernmental Panel on Climate Change (IPCC, 2022) projections suggest that an increase in the rate of Arctic warming is plausible, Bowman Glacier might not survive that long.

Hydrological implications

The reduction in the extent of Bowman Glacier has had a direct impact on the water supply for the Parks Canada basecamp at Tanquary Fiord. Aerial photography and satellite imagery show the presence of glacier ice in the western drainage basin, which drains into May Creek and provides meltwater to the camp, until ~2005 (Figs 3, 7). Since then, streamflow in May Creek has been mainly fed by snowmelt together with active layer and permafrost thaw. Satellite imagery since 2000 shows that in some summers (e.g. 2011, 2015, 2020), both the glacier surface and its three drainage basins become entirely free of snow (Fig. 3). In other years, including at the time of the air photo survey in 2018, seasonal snowpatches persisted throughout the melt season, particularly in narrow depressions at higher elevations where drifting snow tends to accumulate.

In glacierised catchments surface hydrology is modulated by water storage within the glacier system which acts as a buffer, delaying peak seasonal runoff and contributing to streamflow throughout the melt season. Since the loss of glacier ice, the flow regime in May Creek has become dominated by snowmelt, which typically means a shift in the timing of peak discharge, shorter duration of high meltwater input, and increased interannual variability in total runoff. With maximum runoff occurring over a shorter time span earlier in the melt season than before, the rest of the summer is likely to see diminished streamflow and increased influence of rainfall events (Fountain and Tangborn, 1985). May Creek may therefore not remain a reliable source of freshwater for Tanquary Camp in the future. A higher relative contribution of melting ground ice from permafrost degradation is in turn likely to result in more turbid streamflow due to higher sediment, particulate, and solute content (IPCC, 2022).

Conclusion

In this study, SfM-MVS photogrammetry techniques were applied to historical aerial photography from 1959 to create a detailed reconstruction of the surface topography of Bowman Glacier. This enabled direct cloud-to-cloud comparison with recent reconstructions of the glacier generated from aerial photographs acquired in 2018, and calculation of ice surface elevation change and geodetic mass balance spanning six decades. To our knowledge, this study is the first to use historical aerial photography in this way in the Canadian Arctic, allowing for an important gap to be filled concerning changes in glacier surface elevation prior to the satellite record. Although demonstrated here on a single small glacier, the same method can be applied to significantly larger ice masses and on a regional scale.

Over the last six decades, Bowman Glacier lost most of its area ($\sim$78% loss between 1959 and 2020) and volume ($\sim$60% loss between 1959 and 2018). In 2018, the glacier occupied an area of 0.82 km² and had an average ice thickness of 34.1 m. Projecting the rates of change in ice extent and SMB observed over the last decade suggests that the glacier will disappear entirely between 2030 and 2060, or perhaps even sooner. This highlights the importance of ice thickness and volume datasets for documenting past glacier fluctuations and for predicting future changes on a regional basis in the CAA.

Modelling results show the interannual mass balance variability to be dominated by summer balance, meaning that the primary driver for the rapid demise of Bowman Glacier is its sensitivity to increased summer melt due to warmer air temperatures, rather than changes in winter accumulation. Given projections of continued warming across the Arctic, the ongoing trend of rising ELAs is expected to continue, expanding the ablation area of glaciers in the region and leading to increased mass loss. Small ice masses at lower elevations covering narrow elevation spans are particularly sensitive, as even slight changes in ELA can shift entire glaciers into the ablation zone. Those glaciers are clearly out of equilibrium with current climate and, similar to Bowman Glacier, are undergoing committed mass loss and destined to disappear in the future.