Use of historical elevation data to calculate surface-elevation and volume changes of Ha-Iltzuk Icefield, southwest British Columbia, Canada, 1970 to mid-1980s

Abstract

Investigations into glacial changes, including understanding variations in the rates of glacial volume and surface-elevation changes, have increased over the past decade. This study uses historical glacier elevation data in the form of topographic maps from 1970 and a digital elevation model from the mid-1980s to calculate surface-elevation and volume changes for Ha-Iltzuk Icefield, southwest British Columbia, Canada. Results indicate that the icefield thinned at an average rate of 0.76±0.25 ma^–1 during this period. A previous study of Ha-Iltzuk Icefield also using the geodetic method found a thinning rate of 1.0±0.20ma^–1 between the mid-1980s and 1999, indicating a slight increase in the amount of icefield thinning. Within the ablation zone, thinning increased with decreasing elevation at a rate of 1.9±0.68 ma^–1 km^–1 between these two periods (1970 to mid-1980s versus mid-1980s to 1999). Analysis of meteorological data suggests that increases in both temperature and rainfall, as well as decreases in snowfall, likely contributed to the increased thinning rate.

Introduction

Investigations into glacier changes and the relationship of these changes with respect to climate change have greatly increased in the past decade (Dyurgerov and Meier, 1997; Arendt and others, 2002; Rignot and others, 2003; Krabill and others, 2004; Larsen and others, 2007). While much of the concern over sea-level rise is focused on the rapid negative mass imbalances of the Greenland (Zwally and Giovinetto, 2001; Krabill and others, 2004; Luthcke and others, 2006; Rignot and Kanagaratnam, 2006) and Antarctic (Giovinetto and Zwally, 2000; Wingham and others, 2006; Rignot and others, 2008) ice sheets, the current sea-level rise contribution from mountain glaciers may be as significant as the contribution from ice sheets (Kaser and others, 2006). Rignot and Kanagaratnam (2006) estimate the Greenland ice sheet contributed 0.6 mma^–1 to sea-level rise (2000–05), while Kaser and others (2006) estimate all of the world’s mountain glaciers contributed 0.8 mma^–1 to sea-level rise (2001–04). One projection of future mountain glacier and ice-cap changes (excluding the Greenland and Antarctic ice sheets) suggests a loss of 21∓6% of ice volume by the year 2100 (Radić and Hock, 2011).

Despite the large number of studies documenting glacier changes over the past 50 years, only a small percentage of glaciers worldwide have long-term mass-balance records (Dyurgerov and Meier, 1997). The large number of mountain glaciers (~160 000) (Meier and Bahr, 1996) prohibits the possibility of collecting mass-balance data for each glacier. Therefore, many estimates of glacier changes are extrapolated from one representative glacier to all the glaciers in a region (e.g. Fountain and others, 1997). While extrapolation techniques may be necessary to estimate global-scale changes of mountain glaciers, the accuracy of these estimates is limited. Detailed assessments of glacier changes in individual drainage basins are important to validate these large-scale extrapolations.

Remote-sensing techniques are becoming increasingly important in glacier studies, reducing the need for large-scale extrapolations. One important remote-sensing technique is the use of air photos and photogrammetry to create elevation datasets, such as topographic maps and digital elevation models (DEMs). Sets of glacier surface-elevation data collected during different time periods can be used to determine changes in ice surface elevation and ice volume (Arendt and others, 2002; Muskett and others, 2003; VanLooy and others, 2006). This process of mapping glacier surface-elevation changes is known as the geodetic method.

The use of historical data puts current glacier changes in a broader context, allowing us to determine longer-term behavioral trends. Due to their inaccessibility, very few glacierized regions have elevation records that precede the satellite era. In regions where historical air photos exist, calculations can usually be performed on several glaciers in the region, allowing for a more comprehensive assessment of glacier behavior. This study utilizes historical topographic maps created from air photos collected in 1970, as well as DEMs which were created directly from air photos (no intermediate topographic maps created) collected in the mid-1980s, to calculate surface-elevation and volume changes for 30 glaciers which comprise the Ha-Iltzuk Icefield, southwest British Columbia, Canada (Fig. 1). These results are then compared with more recent elevation change results of Ha-Iltzuk Icefield to determine rates of change over time.

Fig. 1. Study area of Ha-Iltzuk Icefield. Numbers show the locations of the 30 Ha-Iltzuk Icefield glaciers which were individually kriged into DEMs from their elevation points (29 numbered glaciers plus Klinaklini Glacier). Bold dates (1984, 1985 and 1986) indicate individual aerial-photo acquisition areas. Dashed lines indicate excluded ice masses which were analyzed by VanLooy and Forster (2008). Solid narrow lines on the icefield indicate locations which were adjusted due to elevation errors within the accumulation area of the DEM.

Study Area

Ha-Iltzuk Icefield is located in the Coast Range of southwest British Columbia, ~300km northwest of Vancouver and ~130km east of the Pacific Ocean (Fig. 1 inset). The 30 glaciers analyzed within this study have an area of 799 km² (based on 1970 topographic maps provided by the Centre for Topographic Information, Canada). This area includes the glaciers draining the icefield, along with eight glaciers adjacent to the icefield. Most of the isolated glaciers and ice masses surrounding the icefield were not included due to problems in delineating their boundaries from the 1970 maps. Some of the isolated glaciers not found in this study were included in a previous study measuring change over a more recent period (VanLooy and Forster, 2008).

Klinaklini Glacier is the largest glacier in the Ha-Iltzuk Icefield; it is ~30 km long and drains nearly half of the icefield area. Several smaller glaciers extend from the icefield, varying in length, with most only a few kilometers but several others reaching 10–12 km long. The mean surface elevation of the icefield in 1970 was 1746 m, with Klinaklini Glacier reaching down to 150ma.s.l., very low for a glacier at 518 N. By 2001 Klinaklini Glacier terminated into a large proglacial lake, but in 1970 this lake had barely begun forming. The other glaciers around Ha-Iltzuk Icefield are land-terminating.

Digital Elevation Model Acquisition and Creation

1970 DEM

Topographic maps of Ha-Iltzuk Icefield provided by the Centre for Topographic Information in Canada were created from air photos taken in 1970 (Natural Resources Canada (NRC), http://maps.nrcan.gc.ca/topo_metadata/topo_metadata_e.php). These maps correspond vertically to the Canadian Geodetic Vertical Datum of 1928 (CGVD28) and horizontally to the North American Datum of 1927 (NAD27), with a contour interval of 100 ft (~30.5 m). The vertical accuracy of each map is 20 m, with a horizontal accuracy of 50 m (NRC, http://maps.nrcan.gc.ca/topo_metadata/topo_metadata_e.php). This horizontal accuracy results in an error in areal extent of ±100m² for the glacier boundaries (Hall and others, 2003). Because the Centre for Topographic Information did not produce corresponding DEMs for these topographic maps, it was necessary to create DEMs for the 1970 glacier elevations to compare with the mid-1980s glacier elevations which were already in DEM format. To produce a DEM of the entire icefield, individual DEMs were created for the 22 glaciers making up the icefield along with the 8 adjacent glaciers (Fig. 1).

DEMs of the 1970 elevations for each of the 30 glaciers were created through the geostatistical interpolation process of kriging the contour point elevations from the topographic maps. Every contour on each glacier was digitized into a series of points (adjusted from feet to meters) spaced no more than 150m apart horizontally. Only the elevation points on the glaciers were used, and no glacier had <348 points. As part of the kriging process, the semivariance of the data points needed to be estimated with a standard statistical model. For consistency, all of the DEMs were created using a Gaussian semivariogram model (which was the best-fitting model for the sample points) with a lag size of 30m and number of lags equaling 25. The only difference in processing each DEM was the search angle for the orientation of the semivariogram model, which is influenced by the flow direction of the glacier. Therefore, the search angle was adjusted so that the number of sampled points in the semivariogram would best fit the Gaussian model for each glacier (Hock and Jensen, 1999; Johnston and others, 2001). The final DEMs were created at a pixel size of 0.7500 (approximately 15m (latitude)±23m (longitude)), with vertical accuracies of the 30 individual glaciers ranging between 20 and 30 m. Finally, the 1970 DEM was horizontally adjusted from NAD27 to NAD83 for consistent comparison with the mid-1980s DEM.

Mid-1980s DEM

A DEM from the mid-1980s was obtained from the Canadian Center for Topographic Information (CTI, http://www.geobase.ca/). The DEM was originally created by the Terrain Resource Information Management (TRIM) program which used air photographs from 1984–86 (dates specific to the study area) to create the DEM directly from photogrammetric interpretation and not from topographic maps (Province of British Columbia, 1992) (Fig. 1). The original digital maps were compiled at a scale of 1 : 20 000 with an absolute horizontal accuracy of ±10m and an absolute vertical accuracy of ±5m (Province of British Columbia, 1992). The DEM was rescaled by TRIM to 1 : 50 000, with the same error estimates, but also containing the same vertical control errors found in topographic maps. The DEM is projected on the NAD83 horizontal coordinates and CGVD28 with a spatial resolution of 0.7500 (same as the 1970 DEM).

Offsets and errors

An analysis of vertical offsets and error between the 1970 and mid-1980s elevations was conducted using six locations chosen off the icefield. DEMs for the offset analysis areas were kriged from the 1970 contour elevation points to compare with the mid-1980s DEM. Bare-earth (unvege-tated) nonglacial areas (as determined by the 1970 topographic maps) with moderate slopes of 10–25˚ (representative of the icefield slopes) were compared. Results of the analysis indicate an offset of 4.4~11 (1σ)m (23 819 total points), with the 1970 DEM having lower elevations. Due to the resulting offset, the 1970 DEM was adjusted by increasing the elevations by 4.4 m. The 1σ (∓11 m) was later used in part to calculate the estimated random error of the difference between the DEMs for the glacier surface-elevation changes.

In flat bright areas such as accumulation areas of icefields, low photographic contrast leads to difficulty in determining accurate elevations, known as ‘contour floating’ (Arendt and others, 2006; Schiefer and others, 2007). The 1970 kriged DEM is likely to contain these errors, as the errors propagate to the DEM because the source data were air photos. The mid-1980s TRIM DEM has also been noted to contain significant errors within the accumulation area (personal communication from M. Milligan 2006). These errors have been described as high-frequency noise in the form of ‘peaks’ and ‘valleys’ running in a north–south direction, as well as ‘bulges’ and ‘dips’ where large sections (tens of km² in area) of the accumulation area rise and drop in a pattern that is not consistent with any realistic environmental variations as compared with other elevation data for the same area during different time periods (VanLooy and Forster, 2008; Fig. 1). To adjust these errors, the high-frequency noise was reduced with a low-pass filter, and the anomalous elevation areas were filled in with the average elevation of the surrounding DEM pixels for the individual time periods (VanLooy and Forster, 2008). Although these errors appeared to occur only in the mid-1980s TRIM DEM, the process of adjustment was conducted for both DEMs for consistency.

Once the DEMs were adjusted, the overall error of surface elevation and volume change was calculated for each of the 30 glaciers. First, random errors were assumed to be the ±1σ of the offset error (11 m) between the two DEMs. The ±1σ was then used to calculate a volume for the entire icefield area. Next, potential volume errors due to the errors in glacier area (±100m²) were taken into account by also calculating a volume using the ±1σ of the offset error. These two volume calculations were then summed in quadrature and the resulting volume was then differenced with the volume change results for the icefield. This produced final volume and, subsequently, surface-elevation change error estimates. The final overall error estimation was ±0.25ma^–1 for the combined 30 glaciers.

Glacier Change Results and Comparison With Other Mountain Glaciers

Calculation of icefield volume change was conducted by differencing the DEMs from the two time periods and multiplying the elevation change of each pixel by the pixel resolution, which were then summed. This volume calculation is described as

(1)

where V is volume change, d ₁ is the 1970 DEM, d ₂ is the mid-1980s DEM, i are the individual pixels for either DEM at location xy, and A is the area of each pixel, i. The surface elevations were then calculated by dividing the volume change results by the maximum extent of the icefield area between the two time periods. Per annum change rates were calculated for each icefield area corresponding with the specific aerial photo acquisition dates for the mid-1980s. Results of the surface elevation calculations indicate that Ha-Iltzuk Icefield changed by an average of –0.76∓0.25 ma^–1 between 1970 and the mid-1980s, which corresponds to a volume change of –9.05±3.01km³ over the ~15 year time period (Fig. 2). The majority of this volume change occurred along Klinaklini Glacier, with a total of – 6.05±1.68km³, which emphasizes the importance of this glacier to the overall mass balance of the icefield (Table 1). Although the overall surface elevation change of the icefield was negative, there were four glaciers that increased in surface elevation and volume. Thickening glaciers are not unexpected, as a previous study has noted positive glacier changes in terms of terminus advances of several glaciers surrounding Mount Waddington, directly to the southeast of Ha-Iltzuk Icefield, during this same period (VanLooy and Forster, 2008).

Table 1. Surface-elevation change rates, ΔZ, of the 30 Ha-Iltzuk Icefield glaciers, and volume change, ΔV, for the period 1970 to mid-1980s; area as determined from 1970 topographic maps; mean surface elevation from 1970 topographic maps; and mean surface elevation from mid-1980s TRIM DEM. Minus values indicate thinning

Fig. 2. Surface-elevation change map of Ha-Iltzuk Icefield, 1970 to mid-1980s. Inset shows a histogram of surface elevation change, ΔZ (positive thickening, negative thinning), by area, A. Solid gray line indicates zero surface elevation change. Dashed gray line indicates mean surface elevation change rate (–0.76 ma^–1).

A previous study of Ha-Iltzuk Icefield compared the mid-1980s TRIM DEM data with the Shuttle Radar Topography Mission (SRTM) DEM data from 2000 (effective 1999) to determine surface-elevation and volume changes between the mid-1980s and 1999 (VanLooy and Forster, 2008). Results from that study show that Ha-Iltzuk Icefield changed at a rate of –1.0±0.20ma^–1. By including the 1970 DEM, we find that the Ha-Iltzuk Icefield experienced a slight increase in thinning rate between the two periods (1970 to mid-1980s: 0.76∓0.25 m a^–1; mid-1980s to 1999: 1.0±0.20ma^–1). However, it should be noted that due to the error estimates this increase in thinning rate may not be significant.

Despite the lack of continuous mass-balance records of glaciers around the world, at least one glacier in southwest British Columbia was continuously measured between the 1960s and the mid-1990s. Place Glacier, located 260 km southeast of Ha-Iltzuk Icefield, experienced a surface elevation change rate of –0.70ma^–1 between 1970 and 1985 (Dyurgerov and Meier, 2004; Fig. 1 inset). However, between 1985 and 1997 the rate nearly doubled to –1.30ma^–1 (Moore and Demuth, 2001). South Cascade Glacier, north central Washington, USA, also has a continuous mass-balance record, which shows that it experienced a more than doubling of surface elevation lowering, from 0.40ma^–1 between 1970 and 1985 to 0.90ma^–1 between 1985 and 1997 (Dyurgerov and Meier, 2004; Fig. 1 inset). While the thinning rate of Ha-Iltzuk Icefield did not double over this same period, it does follow the pattern of increased thinning rates with other glaciers in the region since 1970.

Scatter-Plot Analysis of Surface-Elevation Change Rates

Scatter plots of surface elevation change versus elevation were produced from >140 000 points for the entire icefield and for the two time periods for a more detailed analysis of the surface-elevation change rates (Fig. 3). Preliminary visual analysis indicates that the more recent period (mid-1985 to 1999) had a higher rate of surface elevation change per elevation in the ablation area as compared to the previous period (1970 to mid-1985). For a more in-depth analysis of the surface-elevation changes of the ablation area between the two time periods, it was necessary to estimate the equilibrium-line altitude (ELA) for the icefield so as to separate the icefield into accumulation and ablation areas.

Fig. 3. Scatter plots of surface elevation change, ΔZ, versus elevation, Z, for Ha-Iltzuk Icefield: (a) 1970 to mid-1980s and (b) mid-1980s to 1999. Trend lines exist for points below the estimated ELA, along with respective equations and R ² values. Vertical black line indicates estimated ELA.

The general shape of elevation contours on glaciers changes from concave in the ablation area to convex in the accumulation area, with the inflection contour representing an estimated ELA, which has been shown to be a particularly useful method when using historic elevation data (Leonard and Fountain, 2003). Therefore, contour elevations (100m intervals) were produced for each of the DEMs, and the average icefield ELAs were determined from each glacier’s contour inflection. The estimated average ELA for 1970 was 1484 m, and for the mid-1980s was 1554 m. However, as the analysis of the two scatter plots compares surface-elevation changes over three dates, it is necessary to calculate the average ELA between each date to produce an average elevation boundary between the accumulation and ablation areas. Therefore, the average ELA from the 1999 DEM used in VanLooy and Forster (2008) was also determined, with a result of 1590 m. The average ELA between each period was then calculated to be 1519m between 1970 and the mid-1980s, and 1572 m between the mid-1980s and 1999; these ELA values are used in the scatter plot analysis.

The scatter plot of the accumulation area for both time periods (1970 to mid-1980s, and mid-1980s to 1999) shows thickening and thinning regardless of elevation, with little change overall in either period (Fig. 3). The average surface elevation change in the accumulation area was –0.57±0.16ma^–1 for 1970 to the mid-1980s, but decreased to –0.31±0.20ma^–1 for the mid-1980s to 1999. In the ablation area, thinning rates between the two periods increased substantially. Linear regression trend lines and slope equations of the ablation areas indicate that the rate of surface thinning increased with decreasing elevation during the second period as compared to the first. Between 1970 and the mid-1980s, surface elevation lowered at a rate of 3.2 ±0.71ma^–1km^–1 (Fig. 3a), whereas between the mid-1980s and 1999 it lowered at a rate of 5.1±0.65 ma^–1 km^–1 (Fig. 3b). This represents a substantial increase in sensitivity of ablation rate to elevation.

Comparison With Meteorological Data

To assess the potential causes of the rapid increases in glacier thinning rates since the 1970s for Ha-Iltzuk Icefield, temperature, rainfall and snowfall data were averaged and analyzed for two meteorological stations: Bella Coola to the north (90km from the icefield; 18ma.s.l.) with coastal conditions, and Tatlayoko Lake to the east (94km from the icefield; 870 ma.s.l.) with continental conditions (Adjusted Historical Canadian Climate Data, http://www.cccma.bc.ec.gc.ca/hccd/index.shtml, last accessed June 2010) (Fig. 1 inset). While both meteorological stations are nearly 100 km from the icefield, the data provide a general regional understanding of climate from 1930 to 2000 including the observed time period 1970–2000. Figure 4 (modified from VanLooy and Forster, 2008) shows rainfall, snowfall and temperature between 1930 and 2000 along with the 10 year running mean for each category. Temperatures in the early 1970s are lower than the mean of 5.9˚ C (1930–2000), but gradually rise to >6˚ C by 1985, with an increasing trend into the 1990s. Rainfall appears generally steady throughout the 1970s but begins to increase by the mid-1980s. Snowfall, while high in 1970, declines rapidly by 1985, but rebounds slightly during the 1990s.

Fig. 4. Meteorological graphs for the Ha-Iltzuk Icefield area derived from Bella Coola and Tatlayoko Lake stations: (a) annual rainfall, (b) annual snowfall and (c) annual average temperature. Thick black curves indicate the 10 year running mean. Gray strips indicate observed time period from 1970 to the mid-1980s (effectively 1985). Modified from VanLooy and Forster (2008).

The shift from cooler to warmer temperatures, as well as the shift from higher to lower snowfall, is a possible significant cause of the increased thinning rate of Ha-Iltzuk Icefield. The shift from cooler to warmer temperatures relates to the Pacific Decadal Oscillation, which has been observed in other parts of the Pacific Northwest of North America during the late 1970s (Frauenfeld and others, 2005). Others have observed that increased winter temperatures in this region of British Columbia from 1977 to 1992 have reduced snowpack and increased precipitation as rain instead of snow (Moore and McKendry, 1996). This could explain the overall decrease in annual snowfall as well as the change in thinning rates between the two periods for Ha-Iltzuk Icefield. However, the dominance of the icefield volume change by the 30 km long Klinaklini Glacier could delay the timing of the icefield response to the climatic shift due to the unquantified response time of this glacier.

Summary and Conclusions

Ice surface elevation changes are quantified by comparing a DEM generated from kriging topographic contours with a DEM produced directly from elevations derived through photogrammetry for Ha-Iltzuk Icefield. Results show an average surface elevation change of –0.76±0.25ma^–1 for the period 1970 to the mid-1980s. In comparison, results from a similar study of Ha-Iltzuk Icefield indicate an average surface change of –1.0±0.20ma^–1 for the period mid-1980s to 1999. This increase in thinning rate corresponds to the observed climatic shifts from cooler to warmer temperatures, as well as less snowfall and more rainfall since the early to mid-1980s, but more analysis is needed to quantify the relationship between the climate shift and icefield thinning due to the unquantified response time of Klinaklini Glacier which dominates the icefield volume changes.

A more detailed analysis of icefield thinning was conducted by comparing scatter plots for the accumulation area and ablation areas for the two time periods (1970 to mid-1980s versus mid-1980s to 1999). The average surface elevation change in the accumulation area of Ha-Iltzuk Icefield decreased slightly by 0.26 ma^–1 (from –0.57±0.16 to –0.31±0.20ma^–1). However, the rate of thinning with change in elevation in the ablation area increased from 3.2±0.71 to 5.1 ±0.65ma^–1 km^–1 between the two time periods. The increase of 1.9±0.68ma^–1km^–1 between the two periods pleads for more research on current thinning rates in relation to elevation. This change may be indicative of a dynamic response; however, there are no published velocity data for the Ha-Iltzuk Icefield to measure such a change.