Introduction
Investigations into glacier changes and the relationship of these changes with respect to climate change have greatly increased in the past decade (Reference Dyurgerov and MeierDyurgerov and Meier, 1997; Reference Arendt, Echelmeyer, Harrison, Lingle and ValentineArendt and others, 2002; Reference Rignot, Rivera and CasassaRignot and others, 2003; Reference KrabillKrabill and others, 2004; Reference Larsen, Motyka, Arendt, Echelmeyer and GeisslerLarsen and others, 2007). While much of the concern over sea-level rise is focused on the rapid negative mass imbalances of the Greenland (Reference Zwally and GiovinettoZwally and Giovinetto, 2001; Reference KrabillKrabill and others, 2004; Reference LuthckeLuthcke and others, 2006; Rignot and Kanagaratnam, 2006) and Antarctic (Reference Giovinetto and ZwallyGiovinetto and Zwally, 2000; Reference Wingham, Shepherd, Muir and MarshallWingham and others, 2006; Reference RignotRignot and others, 2008) ice sheets, the current sea-level rise contribution from mountain glaciers may be as significant as the contribution from ice sheets (Reference Kaser, Cogley, Dyurgerov, Meier and OhmuraKaser and others, 2006). Reference Rignot and KanagaratnamRignot and Kanagaratnam (2006) estimate the Greenland ice sheet contributed 0.6 mma–1 to sea-level rise (2000–05), while Reference Kaser, Cogley, Dyurgerov, Meier and OhmuraKaser 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 (Reference Radić and HockRadić 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 (Reference Dyurgerov and MeierDyurgerov and Meier, 1997). The large number of mountain glaciers (~160 000) (Reference Meier and BahrMeier 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. Reference Fountain, Krimmel and TrabantFountain 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 (Reference Arendt, Echelmeyer, Harrison, Lingle and ValentineArendt and others, 2002; Reference Muskett, Lingle, Tangborn and RabusMuskett and others, 2003; Reference VanLooy, Forster and FordVanLooy 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.
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 Reference VanLooy and ForsterVanLooy 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 km2 (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 (Reference VanLooy and ForsterVanLooy 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 ±100m2 for the glacier boundaries (Reference Hall, Bayr, Schöner, Bindschadler and ChienHall 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 (Reference Hock and JensenHock and Jensen, 1999; Reference Johnston, Ver-Hoef, Krivoruchko and LucasJohnston 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 (Reference VictoriaProvince 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 (Reference VictoriaProvince 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’ (Reference ArendtArendt and others, 2006; Reference Schiefer, Menounos and WheateSchiefer 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 km2 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 (Reference VanLooy and ForsterVanLooy 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 (Reference VanLooy and ForsterVanLooy 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 (±100m2) 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
where V is volume change, d 1 is the 1970 DEM, d 2 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.01km3 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.68km3, 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 (Reference VanLooy and ForsterVanLooy and Forster, 2008).
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
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 (Reference VanLooy and ForsterVanLooy 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 (Reference Dyurgerov, Meier, Bamber and PayneDyurgerov and Meier, 2004; Fig. 1 inset). However, between 1985 and 1997 the rate nearly doubled to –1.30ma–1 (Reference Moore and DemuthMoore 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 (Reference Dyurgerov, Meier, Bamber and PayneDyurgerov 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.
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 2 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 (Reference Leonard and FountainLeonard 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 Reference VanLooy and ForsterVanLooy 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 Reference VanLooy and ForsterVanLooy 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.
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 Reference VanLooy and ForsterVanLooy 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 (Reference Frauenfeld, Robert and MannFrauenfeld 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 (Reference Moore and McKendryMoore 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.
Acknowledgements
This work was partially supported by NASA grants NAG5-12552 and NNX09AQ81G, and a Graduate Research Fellowship from the University of Utah. We also thank two anonymous reviewers and the editors for their insightful comments which greatly improved the paper.
