Introduction

Mountain glaciers around the world are rapidly thinning and retreating, and are contributing to sea-level rise (Reference MeierMeier and others, 2007). Most of this retreat is likely due to climatic changes. A recent study of mountain glaciers suggests that more than half the world sea-level rise contribution by meltwater comes from mountain glaciers (∼0.8 mm a⁻¹; Reference Kaser, Cogley, Dyurgerov, Meier and OhmuraKaser and others, 2006). While mountain glaciers are noted to be rapidly changing, many of the calculations for glacier contributions to sea-level rise are extrapolated from calculated ground-based mass-balance changes of nearby glaciers. This extrapolation is often performed due to the difficulty of collecting mass-balance data in situ from a large number of glaciers which are often found in remote locations.

The use of remote-sensing data is an alternative method for collecting data on glacial changes. This can range from optical remote sensing to determine the change in glacial area over time (e.g. Reference Sidjak and WheateSidjak and Wheate, 1999; Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2004), to the use of laser altimetry and digital elevation models (DEMs) created from interferometric synthetic aperture radar (InSAR) as well as from optical Système Probatoire pour l’Observation de la Terre (SPOT) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images to determine surface elevation changes (Reference Arendt, Echelmeyer, Harrison, Lingle and ValentineArendt and others, 2002; Reference Muskett, Lingle, Tangborn and RabusMuskett and others, 2003; Reference AbdalatiAbdalati and others, 2004; Reference Berthier, Arnaud, Vincent and RémyBerthier and others, 2006; Reference VanLooy, Forster and FordVanLooy and others, 2006; Reference Howat, Joughin and ScambosHowat and others, 2007). The major advantage of remote-sensing methods is that it is possible to determine glacial changes over broad regions which are difficult to access. Disadvantages include a coarse temporal resolution and inability to calculate direct mass changes. Also, some remote-sensing techniques, such as those using lidar, still require extrapolation of missing data (Reference ArendtArendt and others, 2006).

This study uses DEMs from the mid-1980s and 1999 to determine the glacial changes of five remote icefields in southwestern British Columbia, Canada (Fig. 1). Shuttle Radar Topography Mission (SRTM) DEMs from February 2000 were compared with Terrain Resource Information Management program (TRIM) DEMs, obtained from the Canadian Center for Topographic Information, derived from stereo air photography flown between 1984 and 1988. A comparison of these two DEMs yields a geodetic mass balance, which in turn can be used to calculate an estimated potential contribution of meltwater to sea-level rise. A recent study by Reference Schiefer, Menounos and WheateSchiefer and others (2007) computes glacier volume change for all glaciers in British Columbia, using SRTM and TRIM DEMs. Our study focuses on a select set of five icefields in the southern coastal region of British Columbia and addresses known errors in the TRIM DEMs. We also include measurement of terminus retreat over three time periods beginning in 1926.

Fig. 1. Study area of southwest British Columbia icefields.

Study Area

The five icefields observed in this study are located in southwest British Columbia, with the most southerly icefield being about 200 km north-northwest of the city of Vancouver (Fig. 1). The first known documented exploration of this remote region was conducted by W.D.A. Munday from 1925 to 1927 (Reference MundayMunday, 1928). This expedition focused on the Mount Waddington glacial area, during which several observations were made of the surrounding glacial landscape and a rough map of glacier termini was compiled. This allows for a comparison of glacier changes with the more recent data, which are discussed later. The four other icefields in this study surround Mount Waddington: the Monarch and Ha-Iltzuk Icefields to the northwest, and the Homathko and Lillooet Icefields to the southeast (Fig. 1). All five icefields have large (10–20 km) valley glaciers extending from their source areas. The glaciers terminate on land (valley-fill deposits) or in proglacial lakes.

Derivation of Digital Elevation Models

To conduct surface-elevation-change and subsequent volume-change calculations, TRIM DEMs were differenced with SRTM DEMs. However, the DEMs were created using different methods, as well as different horizontal and vertical datums, and require an adjustment for consistent comparison. The TRIM DEMs were obtained from the Canadian Center for Topographic Information (http://www.geobase.ca, last accessed February 2008) which were created using air photos from 1984–88 (dates specific to each study area) directly from stereoscopic interpretation (Government of British Columbia, 1992). The original digital maps were compiled at a scale of 1 : 20 000 with an absolute horizontal accuracy of ±10 m and an absolute vertical accuracy of ±5 m (Government of British Columbia, 1992). The DEMs were originally adjusted to North American Datum 1983 (NAD83) horizontal coordinates and Canadian Geodetic Vertical Datum 1928 (CGVD28) with a spatial resolution of 0.75″, which is a posting of approximately 23 m (latitude) and 15 m (longitude) in southwest British Columbia.

The SRTM flown in February 2000 obtained continuous coverage of elevation data from approximately 60° N to 57° S. This was produced from single-pass C-band InSAR data at 5.6 cm wavelengths (Reference FarrFarr and others, 2007). SRTM DEM data are adjusted to the World Geodetic System 1984 (WGS84) horizontal datum and Earth Geopotential Model 1996 (EGM96) geoid heights. For non-United States SRTM data supplied by NASA’s Jet Propulsion Laboratory, the spatial resolution is 3″, which in the study area of southwest British Columbia is a posting of approximately 90 m (latitude) and 60 m (longitude).

Datum Transformations

The DEMs were readjusted to the WGS84 horizontal datum and NAD83 ellipsoidal vertical datum. First, the Canadian DEMs were adjusted from CGVD28 heights to NAD83 ellipsoid heights using the Height Transformation (HT version 2.0) Calculator (http://www.geod.nrcan.gc.ca/apps/gpsh/gpsh_e.php, last accessed February 2008), and the horizontal datums were adjusted from NAD83 to WGS84. Second, for the SRTM vertical elevations the US National Geospatial-Intelligence Agency (NGA) EGM96 Geoid Calculator (http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/intpt.html, last accessed February 2008) was used to first adjust the elevations from EGM96 geoid heights to WGS84 ellipsoid heights. Next, the Horizontal Time Dependent Positioning Calculator (http://www.ngs.noaa.gov/TOOLS/Htdp/Htdp.shtml, last accessed February 2008) was used to determine the offsets between WGS84 and NAD83 ellipsoid heights.

The calculators used to perform the vertical datum transformations only allowed individual points to be calculated. Therefore, a series of point locations were determined for each icefield area spaced approximately 3′ apart (latitudinally and longitudinally). The datum differences were then calculated for each point, and ordinary kriging was used to interpolate the rest of the datum differences between the points. Once these heights were determined, they were added to the SRTM heights adjusted to WGS84 ellipsoid heights, producing SRTM elevations adjusted to the NAD83 ellipsoid. Finally, the Canadian DEMs were resampled to the same size as the SRTM DEMs (90 × 60 m) for consistency in comparing the DEMs.

Icefield Area Determination

The icefield boundaries were created by manually digitizing the glacier outlines from Landsat Thematic Mapper (TM) and Enhanced TM Plus (ETM+) images. These images were obtained in the early autumn, during which time the glacier boundaries are clear of snow, permitting a clearer determination of extent. Landsat ETM+ images from 21–23 September 2000 and 3 October 2001 were used to determine the most recent icefield areas synchronous with the SRTM DEM data. Due to a lack of available Landsat images from the mid-1980s, Landsat TM images from 1990–92 and a multispectral scanner (MSS) image from 1974 were used to determine average distance between the glacier termini to represent the icefield boundaries for the mid-1980s. However, the eastern half of the Lillooet Icefield and the western half of the Monarch Icefield fell outside of the MSS scene, so only the Landsat TM images from the early 1990s were used to approximate the earlier glacier boundaries. An approximated equilibrium-line altitude (ELA) boundary for each icefield was also determined as observed from the September/ October early 1990s and 2000 Landsat images.

DEM Offsets and Errors

Once the DEMs were adjusted to common vertical and horizontal datums, they needed to be corrected for elevation offset relative to each other. Several studies have noted a significant variance in offset in relation to elevation as well as slope (Reference Berthier, Arnaud, Vincent and RémyBerthier and others, 2006; Reference Surazakov and AizenSurazakov and Aizen, 2006; Reference Schiefer, Menounos and WheateSchiefer and others, 2007). To correct for this offset, non-vegetated areas with varying slopes were analyzed surrounding the entire study area. In determining non-vegetated areas it is important to analyze areas in which vegetation is absent from both time periods that the DEMs represent. Therefore, a normalized-difference vegetation index (NDVI) analysis was conducted on the Landsat TM and ETM+ scenes from the 2000/01 and 1990/92 images. This provided images from the two time periods showing healthy vegetation (mostly positive values) and unhealthy or no vegetation (mostly negative values) (Reference KidwellKidwell, 1990). Areas with NDVI values greater than 0.1 were assumed to be healthy vegetation and therefore were masked out. Along with this, areas of snow, ice and water were masked out, therefore leaving values ranging between −0.2 and 0.1 which were assumed as non-vegetated (Reference KidwellKidwell, 1990). Once the non-vegetated areas for both time periods were determined, they were intersected to have one image with areas of no vegetation for both time periods. The no-vegetation image was then intersected with an image of slope values which were calculated by comparison of elevations of adjoining pixels from the DEMs (spatial resolution of 90 × 60 m). Finally, a preliminary subtraction of the TRIM and SRTM DEMs was conducted for the entire study area, and the non-vegetated areas were analyzed to determine the relative offset (9837 points).

Results of the analysis indicated that the offset was relatively stable below 800 m, with an average offset change of 0.3 m (1000 m)⁻¹. Above 800 m the offset value increased by about 11 m per 1000 m of increasing elevation (Fig. 2a). The results also indicated that the offsets are affected by slope, with an increase of 2.7 m(10°)⁻¹ (Fig. 2b). In both cases, the offsets indicated that the SRTM DEM decreases in elevation in relation to the TRIM DEM. Since variations in both elevation and slope influence the offsets, two multiple linear regressions were performed, one for areas below 800 m and another for areas above 800 m:

These equations were then used to correct for the elevation errors in the SRTM DEM.

Fig. 2. Scatter plots of (a) elevation vs preliminary TRIM – SRTM difference image offsets (trend line equation above 800 m is indicated); and (b) slope vs preliminary TRIM – SRTM difference image offset (trend-line equation below 800 m is indicated).

In the event of horizontal offsets, a dramatic systematic variation in elevation is noticeable in the preliminary difference image along sharp topographic gradients such as mountain peaks and glacier-valley wall boundaries, indicating a horizontal misalignment between the DEMs. However, no systematic shift of this kind was found on the preliminary difference image, so no horizontal offset corrections were performed. This was also the case for Reference Schiefer, Menounos and WheateSchiefer and others (2007).

During the process of creating the TRIM DEMs from the air photographs, systematic errors occur primarily in areas of flat bright terrain, such as in the accumulation areas of the icefields, where it is difficult to determine accurate elevations due to low photographic contrast (personal communication from M. Milligan, 2006). This appears to be evident in the TRIM DEMs, as there are several locations in the accumulation areas where there are strangely varying elevations (as much as 100–200 m over 0.5 km) across a region that should have a much lower slope (∼30 m over 0.5 km) as indicated by analyzing topographic maps (from the 1970s) of these same areas. One method of correcting these errors would be to use a topographic map from the same time period to manually adjust the incorrect elevations (Reference HeipkeHeipke, 1995); however, no paper maps are available for the TRIM time period. While the topographic maps are from time periods of approximately 10 years before the air photos from which the DEMs were created, it is highly unlikely that these large variations in elevation would form over this period. These varying elevation patterns come in two forms (Fig. 3). First, east–west profiles set across the accumulation areas show a regular pattern of high-frequency noise, with elevation changes as much as 60–100 m spaced every 500–700 m. Along with this, there are much more abrupt and spatially extensive elevation variations in formations of ‘bulges’ and ‘dips’ on the order of ±180 m, again in areas that are shown to be flat on the 1970s topographic maps. In one instance on the Monarch Icefield, a bulge of tens of metres occurs directly next to a ‘dip’ also on the order of tens of metres (Fig. 3). These elevation variations happen to occur almost directly south and north of each other as if along a seam of two air photos, with a drastic change in elevation between them. Many other bulges and dips occur across all the icefields, but this case appears to be the worst. The method for addressing these two problems is discussed later.

Fig. 3. TRIM – SRTM elevation differences of the Monarch Icefield, highlighting the high-frequency noise as well as the bulge and dip errors located in the accumulation areas.

In comparing surface elevation changes, it is important to compare data from the two time periods as closely as possible to the same seasons to obtain accurate mass-balance changes. The TRIM DEMs were created from air photos collected in the mid-1980s mostly during the months of July and August, while the SRTM was flown in February 2000. However, the C-band signal from SRTM penetrates through dry snow during the cold winter months such as February, theoretically by as much as >10 m (Reference Ulaby, Moore and FungUlaby and others, 1986), and through observed penetration by as much as 9 ± 2 m (Reference Rignot, Echelmeyer and KrabillRignot and others, 2001). Another study of glaciers in the French Alps shows that C-band SRTM closely matches the September 1999 elevations, particularly in the ablation area, but above 2200 m the SRTM elevations were underestimated by 4.1 ± 1.6 m, most likely due to radar penetration of the snow surface (Reference Berthier, Arnaud, Vincent and RémyBerthier and others, 2006). Therefore, as it is difficult to determine and correct for the depth of radar penetration, it should at least be stated that accumulation-area elevations could have errors of as much as 10 m or more, depending on the density and water content of the snow (Reference Ulaby, Moore and FungUlaby and others, 1986).

The results of the preliminary DEM subtractions of each icefield highlighted the TRIM elevation variation problems in the accumulation areas. This was especially true for the larger bulge and dip variations. While the TRIM DEMs contained these variations, the SRTM DEMs over the same areas were fairly flat, much like the topographic maps of the 1970s. Statistical analysis of the DEM subtractions in the accumulation areas showed that isolating >±1σ of the DEM subtraction nearly completely removed the large variations. However, these areas now had to be adjusted to represent the missing areas with more realistic elevations. Other studies suggest that correcting for errors, such as sinks (isolated pixels) which do not exist in this dataset, can be done by assigning the sink pixels the same elevation as the lowest surrounding pixel elevation (Reference Ledwith and LundénLedwith and Lunden, 2001). However, the area of elevation anomalies in this study can be as extensive as 900+ pixels (approximately 5 km²). Therefore, the average elevation of the surrounding points for each missing area was calculated and used to fill in the missing space. The average elevation for the SRTM DEMs of the points surrounding the missing areas were also calculated and used to fill in the missing area for the SRTM DEM, which applied uniform processing steps to both the TRIM and SRTM datasets.

Next, to reduce the high-frequency noise within the accumulation area, a 9 × 9 low-pass filter was used to smooth the accumulation area for both the TRIM and SRTM DEMs. This was done for visual reasons since a low-pass filter smooths the appearance of the icefield and does not affect the mean volume or surface-elevation change calculations. The systematic problems were only corrected for in the accumulation area, as the ablation area does not contain these problems. While the high-frequency noise does not necessarily occur across all accumulation areas, the entire accumulation area of each icefield was filtered, to be certain to filter all of the problem areas including those difficult to identify.

Once the DEMs were adjusted for systematic problems, a calculation of errors for the surface-elevation and volume changes was conducted. The random errors in the DEMs are assumed to be the ±1σ of the offset between the TRIM and SRTM DEMs, with each 100 m elevation contour as an independent sample point. Thus, the random errors for the surface-elevation changes between the two DEMs were calculated by dividing the ±1σ by the number of 100 m contour intervals for each individual icefield, as used by Reference Rignot, Rivera and CasassaRignot and others (2003) and Reference VanLooy, Forster and FordVanLooy and others (2006). Finally, the error bounds for the systematic errors remaining in the DEMs (after the systematic error corrections) were determined by adding 5 m (absolute TRIM vertical accuracy) to the TRIM DEM for the accumulation area. A new volume change was calculated, and the difference between the new error-estimated volume change and the original volume change was taken as the systematic error. Finally, the sum of the random and systematic errors was taken as the overall error for the volume and surface-elevation changes (Table 1).

Table 1. Average icefield area determined by Landsat images (Ave. A), volume changes (ΔV), area-averaged surface elevation changes (ΔZ), average area below ELA (Ab), volume change below ELA (ΔVb), change in area between the mid-1980s and 1999 (ΔAm), percent of area above ELA identified as anomalous and replaced (% Aanom), mean surface elevation (Mean Z) of the five southwest British Columbian icefields, approximate ELA (Appro. ELA), volume changes without adjustments in the accumulation area (ΔV w.a.), and surface elevation changes without adjustments in the accumulation area (ΔZ w.a.). Minus signs indicate thinning or volume loss

Terminus and Area Changes

The icefield areas determined from the Landsat images allowed for an analysis of the change in surface area between the mid-1980s and 2000/01. These calculations were determined by analyzing the change in the periphery of the icefields, so this analysis does not account for any change in area that may have occurred within the icefields above the ELA. The reason for only examining change in area by analyzing the terminus positions was the difficulty in determining icefield area in the accumulation area due to snow cover on nunataks. The results show that the total icefield area decreased by 51.3 km², with the greatest change occurring on the Ha-Iltzuk Icefield (13.8 km²), mostly on Klinaklini Glacier (Table 1).

Terminus-position changes of some of the larger glaciers were also determined from the MSS and Landsat imagery between 1974 and 1990–92, and between 1990–92 and 2000/01. This was done by taking an average distance of several (six to eight) measurements across the front of the glacier, as a single measurement would not provide a good estimate of terminus position change due to uneven movement of the terminus. Error estimations for the position changes were determined to be equal to one pixel of the images (57 m for the 1974 MSS image, and 29 m for the Landsat TM and ETM+ images).

Fifteen glacier terminus changes, sampled from each of the icefields, were measured, showing there was an average terminus retreat of 290.0 ± 6.0 m (20 ± 3.0 m a⁻¹) from 1974 to 1990–92 (Table 2). Over the period 1990–92 to 2000/01, there was an average (over 19 glaciers) terminus retreat of 430.0 ± 3 m (40.0 ± 3.0 m a⁻¹),with two glaciers retreating more than 1 km. The terminus of Jubilee Glacier was heavily debris-covered in the 1974 MSS image, so a calculation of the terminus change was determined only for the period 1990–92 to 2000/01. Three other glaciers (Fyles, Lillooet and Bridge) were outside the 1974 image and also only had terminus-change determinations for 1990–92 to 2000/01. One particularly large glacier (Tiedemann) was not included in this analysis since the terminus is very heavily covered in debris and no distinct position could be noted on the images. However, for the overall determination of area for the Mount Waddington glaciers, an estimated terminus was used.

Table 2. Sampled glacier terminus changes and retreat rate per year from all five icefields for the periods 1927–74, 1974 to 1990–92 and 1990–92 to 2000/01 for the Mount Waddington area. NA: no data available

With the historical map of the Mount Waddington area glaciers created during the Munday expeditions of 1926 (to the east of Mount Waddington) and 1927 (to the west of Mount Waddington) (Reference MundayMunday, 1928), it was possible to roughly determine the amount of terminus change for several glaciers since the late 1920s (Fig. 4). The map does contain geographic coordinates and several landmarks, such as mountain peaks, rivers and creeks, and even in one case a landslide which is still visible on the Landsat images. These landmarks allowed for the Munday map to be georegistered to the Landsat images and then to have the 1926/27 terminus positions outlined. However, despite the usefulness of the landmarks, the Munday map is still quite rough (terminus changes could be off by as much as ±300 m as determined by averaging the differences of distances between coincident points (e.g. stream intersections) on the Munday map and the Landsat images) in the geographic positioning of the landmarks, as well as the glacier termini, but, as many of the glaciers have retreated several km up-valley, it still provides a useful historical description of the terminus changes. The average terminus retreat between 1926/27 and 1974 was 2490.0 ± 260 m (52.0 ± 6.0 m a⁻¹), with the greatest terminus retreat occurring in Franklin Glacier as it retreated 4100 m (Table 2). Nine other measured glacier terminus changes showed retreat ranging from 1100 to 3200 m. Two glaciers (Jubilee and Confederation) were tributaries of Franklin Glacier in 1927 but have since retreated several km from the confluence. Overall, the average change in terminus positions appears to fluctuate, with both retreats and advances occurring over the entire time period 1926–2001.

Fig. 4. Munday expedition map of Mystery Mountain (Mount Waddington), showing glacier locations and terminus positions from the 1920s. Terminus positions determined from a 2001 Landsat image and glacier names are shown in bold. (Modified from Reference MundayMunday (1928) with permission from Geographical Review.)

Ice-Volume and Surface-Elevation Changes

The total glacial volume change for the five icefields in this study from the mid-1980s to 1999 is −19.4 ± 8.8 km³ (−1.5 ± 0.7 km³ a⁻¹)over an average glacial area (the average area between the mid-1980s and 1999) of 3218 km². This equates to an area-averaged surface elevation change of −6.0 ± 2.7 m (−0.5 ± 0.2 m a⁻¹) (Table 1). Individually the total volume change of the icefields ranged from −11.1 ± 1.8 km³ (Ha-Iltzuk Icefield) to −1.2 ± 1.9 km³ (Lillooet Icefield), with the area-averaged surface elevation changes ranging from −1.0 ± 0.2 m a⁻¹ (Ha-Iltzuk Icefield) to −0.2 ± 0.3 m a⁻¹ (Lillooet Icefield) (Fig. 5). While it is difficult to accurately determine the amount of meltwater equivalent from these icefield volume changes due to lack of knowledge about the actual glacial mass densities across the icefields, most of the volume loss is at the ice fronts, and below the equilibrium line, so it is appropriate to use an ice density of 0.9 g cm⁻³. This results in a meltwater equivalence of −17.5 ± 8.0 Gt (−1.3 ± 0.6 Gt a⁻¹) which in turn equates to an approximate sea-level rise of 0.004 ± 0.002 mm a⁻¹ . A second analysis of total volume change, neglecting the extra correction for variable elevation and slope offsets, used a constant mean offset of 2.6 m between TRIM and SRTM. This increased the estimate of volume loss considerably to 56.0 ± 14.6 km³. However, we believe that the variable-offset approach resulted in more trustworthy numbers.

Fig. 5. TRIM (mid-1980s) and SRTM (effective 1999) DEM difference images for (a) Monarch and Ha-litzuk Icefields; (b) Mount Waddington Glacial Area and Whitemantle Icefield; and (c) Homathko and Lillooet Icefields (inset shows icefield elevation vs area-averaged surface elevation change for each icefield). White lines on the icefields represent the approximate ELA.

Observations of the DEM difference images show a clear pattern of rapid thinning around the periphery of the icefields, with some thickening at higher elevations (Fig. 5). Ha-Iltzuk Icefield, which is the largest icefield and has had the greatest volume change, is particularly interesting, as nearly half of the icefield is drained by one glacier (Klinaklini). This glacier extends down to 165 m a.s.l., which is extremely low for a temperate glacier at approximately 51° N. At the lower elevations, the glacier has thinned by >100 m between the mid-1980s and 1999. Overall, Ha-Iltzuk Icefield has a mean elevation of 1675 m (dominated by Klinaklini Glacier) which is the lowest mean elevation of the five icefields (Table 1). Another glacier to note is Bridge Glacier on Lillooet Icefield (Fig. 5c) which, unlike most of the other outlet glaciers on the icefields, has thickened at lower elevations. This raises questions about the flow dynamics of the glacier, such as surging; however, no published accounts of surging in this area could be found.

To compare surface-elevation change with icefield elevation, the average elevation change was calculated for every binned 100 m of elevation (inset of Fig. 5c). The thickening in the accumulation area of each icefield is apparent mostly above 2000 m. The one exception is Monarch Icefield which shows a general surface lowering in the accumulation area; however, this is likely due to the large errors from the TRIM DEMs, as this icefield contained the largest and most numerous problematic areas in the accumulation zone. Below the ELA, the thinning rates vary between each icefield but generally increase with decreasing icefield elevation. However, at lower elevations the thinning rates vary much more between icefields as compared with the accumulation areas.

Finally, ice-volume and surface-elevation changes were calculated for the five icefields without adjusting for the errors in the accumulation zone. The purpose was to analyze the significance of these errors. For all five icefields, the difference between the adjusted and unadjusted calculations fell within the error bounds, but the unadjusted calculations for all icefields indicated less thinning (Table 1). The total volume change for the unadjusted calculations was −14.9 ± 8.1 km³ vs −19.4 ± 8.8 km³ for the adjusted calculations. This implies the net positive and negative adjustments to the surface are within the error bounds of the volume calculation.

Comparison of Glacial Changes with other North American Icefields and Glaciers

While using remote-sensing data to compare thinning rates and volume changes is informative for understanding glacial changes over large regions, unlike with continuous field-based measurements one major problem that occurs in comparing DEMs is the lack of consistent periods of glacial changes between two different locations. This can cause the comparison of results to be misleading; however, due to the lack of data, there may be no better way than to compare volume and thinning rates for time periods that are as close as possible.

Thinning rates of glaciers and icefields for other parts of British Columbia and Alaska appear to be faster than for the five southwest British Columbian icefields. One comprehensive study of Alaskan glacier melt rates suggested average thinning rates of 0.52 m a⁻¹ between 1950 and the mid-1990s, whereas the rates for those same glaciers increased to 1.8 m a⁻¹ between the mid-1990s and 2000/01 (Reference Arendt, Echelmeyer, Harrison, Lingle and ValentineArendt and others, 2002). The icefields of the southeastern Alaskan panhandle (e.g. Juneau Icefield, Stikine Icefield, Glacier Bay area) (14 570 km²) have been shown to be thinning at an average of 1.0 ± 0.3 m a⁻¹ over the period 1948–79 to 2000 (Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others, 2007). The total volume change for Kenai Peninsula icefields in south-central Alaska between 1950 and 1999 indicates a thinning rate of 0.6 ± 0.1 m a⁻¹ (Reference VanLooy, Forster and FordVanLooy and others, 2006). However, it is important to note again that the time periods for the Kenai Peninsula and southeastern Alaskan panhandle icefields extend further into the past and therefore are likely to be influenced by different climatic conditions which would affect averaged thinning rates. At similar latitudes to those of the five British Columbian icefields, Place Glacier (Fig. 1), located between Lillooet Icefield and Vancouver, thinned at a rate of 1.3 m a⁻¹ between 1985 and 1997, and by as much as 1.4 m a⁻¹ throughout the 1990s (Reference Moore and DemuthMoore and Demuth, 2001), but Peyto Glacier in the Canadian Rocky Mountains has thinned by 0.6 m a⁻¹ (Reference Dyurgerov, Meier, Bamber and PayneDyurgerov and Meier, 2004), also between 1985 and 1997.

A direct comparison with Reference Schiefer, Menounos and WheateSchiefer and others (2007) is not possible, for two reasons. While they subdivide British Columbia into mountain regions, the area containing the five icefields in our study (South Coastal Mountains) includes significantly more glacier area. Their South Coastal Mountains region glacial area is 8037 km², compared with 3218 km² for the five icefields in this study. There is also a difference in the variable offset used to adjust the SRTM DEM. Reference Schiefer, Menounos and WheateSchiefer and others (2007) found an average SRTM offset of −12 m km⁻¹ compared to the −11 m km⁻¹ for this study. In addition, this study variably adjusts the SRTM offset relative to slope. Reference Schiefer, Menounos and WheateSchiefer and others (2007) find a mean thinning rate of −0.9 ± 0.2 m a⁻¹ over the South Coastal Mountain region compared to −0.5 ± 0.20 m a⁻¹ for the five icefields in this study. However, the five icefields in this study have large accumulation areas (relative to the ablation areas), therefore reducing the amount of area-averaged elevation change compared to smaller glaciers in the rest of the region (such as Place Glacier) which might have dominating ablation areas.

Meteorological Data Analysis

An average of the annual temperature anomalies from the 71 year mean (1930–2000), as well as average snowfall and rainfall were analyzed over the same time period from the two closest stations surrounding the icefields (locations in Fig. 1 inset). The data contain few missing observations, but, in the case of missing temperature data, values were estimated from highly correlated neighboring stations (Adjusted Historical Canadian Climate Data, http://www.cccma.bc.ec.gc.ca/hccd/index.shtml, last accessed February 2008). However, for annual snowfall and rainfall, estimations were not made for stations with two or more consecutive missing years.

The 10 year running mean was calculated for all three variables for a comparison of change over time. In general, increasing annual rainfall and decreasing snowfall coincided with increasing temperature. In particular, the 10 year running mean for the annual temperature anomalies from the 71 year mean indicates a strong shift from negative to positive values around 1986 (Fig. 6). The annual snowfall also shows what appears to be a trend to decreasing amounts since the late 1970s, while rainfall indicates an increase during this same period. It is likely that the increase in temperatures along with the decrease in snowfall and increase in rainfall between the mid-1980s and 1999 has contributed to the thinning of the icefields. Reference Moore and McKendryMoore and McKendry (1996) also noticed this change in precipitation for southwestern British Columbia between 1977 and 1992, and particularly the decline of winter snowpack due to increased rainfall during the winter. The averaged terminus retreat rates do not appear to respond in a simple way to temperature anomalies and precipitation trends recorded at the meteorological stations. This is evident as the glaciers retreated most rapidly during the period 1927–74 when temperatures and rainfall were generally lower and snowfall amounts were higher on a 10 year running mean. The opposite was the case for the period 1974–1990/91 which experienced increasing temperatures and rainfall with decreasing snowfall, during which time the glacier terminus retreat slowed.

Fig. 6. Annual rainfall (a), annual snowfall (b) and annual temperature anomalies against the 71 year mean (c), as averaged from two climatological stations. Thick black line indicates the 10 year running mean.

Summary

The five southwest British Columbia icefields discussed in this paper have thinned between the mid-1980s and 1999 by an average of −0.5 ± 0.2 m a⁻¹. Although the total area of the icefields is small compared to some of the Alaskan and northwest British Columbian icefields, the rapid thinning rate has led to a potential meltwater contribution of 1.3 ± 0.6 Gt a⁻¹. In comparison, the five southwest British Columbia icefields are contributing as much as 0.5% of the total world mountain glaciers to sea-level rise between the mid-1980s and 1999 (Reference Kaser, Cogley, Dyurgerov, Meier and OhmuraKaser and others, 2006).

The historical map provided by the Reference MundayMunday (1928) expedition allowed for a rough comparison with 1974 MSS, and 1990s/2000s Landsat images of glacier termini changes for the Mount Waddington area. This comparison indicated that the average speed of terminus retreat was three times faster between 1927 and 1974 (52 ± 6 m a⁻¹) than between 1974 and 1990/91 (17 ± 3 ma⁻¹). Terminus retreat rates then increased again between 1990/91 and 2000/01 (44 ± 3 m a⁻¹).

Results for these icefields indicate care should be taken when accounting for SRTM offsets between other DEMs, particularly at elevations above 800 m. Comparison of volume-change results between SRTM DEMs adjusted with a constant offset and SRTM DEMs adjusted with variable offsets in relation to elevation and slope indicates a potential overestimation of volume change by more than a factor of two when using the constant offset. This emphasizes the importance of adjusting for the variable offsets to provide more accurate measurements of mountain glacier ice-volume change when using SRTM-derived DEMs.