Introduction
Studies of historic glacier change commonly are based on comparison of a digital elevation model (DEM) with an older DEM or topographic map (Reference D’Agata, Smiraglia, Zanutta and ManciniD’Agata and others, 2005; Reference Granshaw and FountainGranshaw and Fountain, 2006; Reference Surazakov and AizenSurazakov and Aizen, 2006; Reference D’Agata and ZanuttaD’Agata and Zanutta, 2007; Reference Jackson and FountainJackson and Fountain, 2007; Reference Racoviteanu, Manley, Arnaud and WilliamsRacoviteanu and others, 2007; Reference Schiefer, Menounos and WheateSchiefer and others, 2007; Reference Wang, Hou, Hong, Hur and LiuWang and others, 2008). The accuracy of the derived elevation and volume changes depends on how accurately the glacier surface is represented in these topographic datasets.
Low photogrammetric contrast in the accumulation areas of alpine glaciers can introduce large errors in DEMs, even using modern software and workstations (Reference Cox and MarchCox and March, 2004). Even recent, non-photogrammetric, digital elevation datasets that allow a more rigorous analysis of error can be biased over glaciers. Reference Schiefer, Menounos and WheateSchiefer and others (2007), for example, showed that DEMs produced for some British Columbia (Canada) glaciers from single-pass interferometric synthetic aperture radar (InSAR) provided by the Shuttle Radar Topography Mission (SRTM) have elevation errors of up to 12 m km−1, approximately twice the error reported for SRTM datasets by others (Reference Berthier, Arnaud, Vincent and RémyBerthier and others, 2006; Reference PaulPaul, 2008).
Given these errors, independent audits of glacier surface elevations are needed to assess the reliability of published maps and DEMs. The European Remote-sensing Satellite (ERS)-1/-2 tandem repeat-pass InSAR dataset, with its near-global coverage in 1995/96, is a valuable resource for producing both DEMs and maps of Earth surface displacement (Reference Goldstein, Engelhardt, Kamb and FrolichGoldstein and others, 1993; Reference Rabus, Lang and Sawaya-LacosteRabus and Lang, 2000), especially in high-latitude areas beyond the range of the 2001 SRTM. Caution is required, however, in using these datasets, because differential, repeat-pass interferometric (D-InSAR) processing techniques (e.g. Reference Joughin, Kwok and FahnestockJoughin and others, 1998) may introduce substantial systematic errors over valley glaciers where topography is unknown and surface velocity changes seasonally.
We describe a new D-InSAR method for constructing an accurate DEM of a glacier surface and apply the method to Black Rapids Glacier, Alaska, USA. We use the new DEM to assess the accuracy of other topographic datasets and to show how Black Rapids Glacier has adjusted since its last surge in 1936/37.
Study Area
Black Rapids Glacier is a 40 km long valley glacier located in the central Alaska Range of interior Alaska (Fig. 1). It has been extensively studied by scientists working for the US Geological Survey, University of Alaska Fairbanks and University of Washington since the 1960s (Reference PostPost, 1960; Reference Harrison, Mayo, Trabant, Weller and BowlingHarrison and others, 1975; Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others, 1996; Reference Truffer, Motyka, Harrison, Echelmeyer, Fisk and TulaczykTruffer and others, 1999; Reference Rabus and FatlandRabus and Fatland, 2000; Reference Amundson, Truffer and LüthiAmundson and others, 2006).
Fig. 1. Map of Black Rapids Glacier showing the trace of the Denali Fault (dashed curve) and distances from the head of the glacier in 5 km increments (modified from Reference Amundson, Truffer and LüthiAmundson and others, 2006). The location of the Susitna–Black Rapids glacier divide is shown.
Black Rapids Glacier is a surge-type glacier. It last surged about 6.4 km in 1936/37, nearly blocking Delta River in the trunk valley below the glacier (Reference HanceHance, 1937; Reference PéwéPéwé, 1951). There is geologic evidence for several late-Holocene surges of Black Rapids Glacier, at least two of which completely blocked Delta River (Reference Reger, Sturmann, Beget, Solie and TannianReger and others, 1993; Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others, 1996). The surge periodicity of Black Rapids Glacier is not known with certainty, although Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996) estimate it to be about 50–75 years based on the size and position of surge loop moraines that are being pushed out by the Loket tributary. It is also unclear whether the pronounced warming of climate that is occurring today in the Alaska Range will terminate the surge cycle of Black Rapids Glacier, as has been the case in other regions (Reference HoinkesHoinkes, 1969; Reference Dowdeswell, Hodgkins, Nuttall, Hagen and HamiltonDowdeswell and others, 1995).
The ablation zone of Black Rapids Glacier lies in an east-trending valley that also contains the trace of the Denali Fault. Most of the accumulation zone is in two north-facing valleys (Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others, 1996), one of which contains the Loket tributary. The Loket tributary supplies substantial ice to Black Rapids Glacier between surges.
Methodology
We constructed a DEM of the 1995 surface of Black Rapids Glacier using ERS-1/-2 repeat-pass satellite radar interferometry. Interferometric phase, φ, is computed as a difference between two SAR acquisitions. It comprises independent contributions from orbit geometry, φ o, surface topography, φ t, surface motion, φ m, atmospheric changes, φ a, and sensor noise, φ n (e.g. Reference Rabus and FatlandRabus and Fatland, 2000):
The topographic contribution is proportional to the baseline, B p, which is the difference in satellite positions between the SAR acquisitions measured perpendicular to the satellite line-of-sight. The motion contribution is proportional to the scalar product of the surface displacement vector d and the satellite line-of-sight vector nLOS (Reference Rabus and FatlandRabus and Fatland, 2000):
The goal of D-InSAR is to isolate φ m by removing all other contributions. Orbital parameters are provided with the SAR data, so φ o can be modeled. φ t can be calculated from an accurate external DEM if it is available; otherwise φ t must be determined simultaneously with φ m, using two or more interferograms with different B p. A variant of this approach is used here. φ a and φ n are generally sources of error that limit the accuracy of φ m and φ t. The latter can be reduced by spatial low-pass filtering, but at the cost of spatial resolution. In some cases, as in this study, φ a can be estimated and removed from the data.
Black Rapids Glacier is north of 60° N, so no SRTM data are available as an external DEM. A high-resolution (10 m) InSAR DEM was produced by Intermap Technologies in 2001, but it only covers the lowermost parts of the ablation area. An old US Geological Survey (USGS) topographic dataset is available for Black Rapids Glacier (USGS, 1951a,b). It was published as a DEM in 1999 but is based on topographic maps (Mount Hayes B-4, B-5 and B-6) published in 1951. Map sheet B-5 covers most of Black Rapids Glacier, and most of the mapping photography was acquired in 1949 (Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others, 1996). We therefore assume that the USGS DEM represents the glacier surface in 1949. The DEM has approximately 60 m resolution and 7–15 m vertical accuracy. Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996) pointed out that the USGS maps may have significant biases, especially over the accumulation area of Black Rapids Glacier, which was a featureless snowfield at the time of photography. They were unable, however, to make a quantitative assessment of the error.
Figure 2 shows the approach we used to analyze the ERS tandem data for Black Rapids Glacier. We created small and large baseline interferograms to (1) determine the degree to which velocity is temporally constant and (2) isolate the phase due to elevation difference with respect to the USGS DEM. After unwrapping this difference phase, we (3) corrected unwrapping errors with a new semi-automatic multi-baseline approach and (4) reduced the atmospheric phase error from all original interferograms through a combination of filtering and interpolation. Lastly, (5) we produced a new DEM in map coordinates as the sum of the USGS DEM and the D-InSAR-derived difference DEM.
Fig. 2. Flow chart for DEM construction. Labels 1–5 refer to steps described in the text.
Figure 3 shows the two differential interferograms involved in step 1 of the process, one with a large baseline (468 m) and the other with a small baseline (8 m). The differential interferograms were produced from three 1 day repeat interferograms and flattened with the USGS DEM. The baselines of the differential interferograms are the differences of the baselines of the original tandem interferograms. For example, the tandem interferograms from November 1995 (12 and 13 November; B p = 166 m) and December 1995 (17 and 18 December; B p = 158 m) yield a differential baseline of 8 m.
Fig. 3. Differential interferograms of Black Rapids Glacier and its surrounding area: (a) topography change only; (b) velocity change only. Black areas are areas of layover, shadow and poor coherence. Note the high density of fringes on the glacier in (a) and the lack of fringes on the glacier in (b).
Differential interferograms contain signals of both temporal velocity change and topographic change, the latter relative to the external DEM. If, however, velocity is constant between the repeat interferograms, the differential interferogram records only topographic change. The small baseline differential interferogram has a negligible topographic signal and can therefore be used to assess the magnitude of velocity change over the acquisition period of the complete tandem dataset. Inspection of the small-baseline differential interferogram shows that velocity changed little between mid-November and mid-December 1995 (Fig. 3b): significantly less than one fringe everywhere on the glacier surface. We conservatively estimate the corresponding velocity change to be <2 cm d−1 along the line-of-sight. This finding is consistent with independent observations of small winter velocity variations of Black Rapids Glacier (Reference Rabus and FatlandRabus and Fatland, 2000).
A consequence of this finding is that we can assume the large-baseline (−468 m) differential interferogram between October 1995 (8 and 9 October; B p =−302 m) and November 1995 is dominated by topography. In the ablation zone of Black Rapids Glacier, where significant surface change is expected with respect to the USGS DEM, the large-baseline differential interferogram (Fig. 3a) shows a high fringe density, whereas in areas where the USGS DEM should be accurate (e.g. in the Delta River valley) the fringe density is much lower.
The next step (2 in Fig. 2) is to unwrap the large-baseline differential interferogram. We improved the coherence with an adaptive spectral filter (Reference Goldstein and WernerGoldstein and Werner, 1998) and then used Reference Costantini, Guyenne and DanesyCostantini’s (1996, Reference Costantini1998) minimum-cost flow algorithm to unwrap the phase using the coherence as costs. The adaptive spectral filter is nonlinear but on average reduces the spatial resolution of the data from 20 m to ∼30 m. Since radar phase values are periodic functions of 2π, they automatically wrap after reaching 2π, for example as a result of surface motion. Phase unwrapping allows retrieval of the phase information, but requires control points of known velocity. The minimum-cost flow unwrap-per minimizes the weighted deviations of unwrapped phase between neighbouring pixels.
The next step (3) is to correct for phase-unwrapping errors originating from fringe aliasing due to strong velocity gradients and numerous noisy (incoherent) areas. The latter result from SAR layover on steep slopes surrounding the glacier. SAR layover occurs where slopes facing the satellite are steeper than the sensor’s look angle. Phase-unwrapping errors are typically manifested as polygons that are an integer multiple, N, of 2π, above or below the background phase surface. If there are numerous error polygons, it may be impossible to recognize and correct them automatically, and it can be difficult to correct them interactively.
Interactive delineation and correction of unwrapping errors depend on the colour range specified by the user to display the interferograms. If the colour range spans the full spectrum of the unwrapped phase values, the contrast may be insufficient to detect and delineate error polygons for which N is a small value, because those error polygons will have about the same colour as the background.
An alternative way of visualizing the error polygons is to remove the topography, which contains the unwrapping errors of the large-baseline interferogram, from an original, wrapped, tandem interferogram. These errors are then revealed as sharply delineated polygons differing from the background colour in the now-flattened, wrapped tandem interferogram. This procedure allows us to manually digitize the error polygons in ArcGIS and to eliminate them by iteratively raising or lowering the polygons in the unwrapped, large-baseline differential interferogram. The correct multiple, N, of 2π is found when the colour differences disappear.
A phase-unwrapping error polygon is only evident in a tandem interferogram if N is different from an integer multiple of the baseline ratio between the large-baseline differential interferogram and the tandem interferogram. For example, the baseline ratio of the large-baseline (B p = −468 m) interferogram to the November 1995 (B p = 166 m) tandem interferogram is ∼3, which means that erroneous polygons of N = ±1 and N = ±2 will be revealed, but polygons of N = ±3 will not be visible. However, if a second tandem interferogram (e.g. 31 March and 1 April; B p = 65 m) is considered, the baseline ratio is about 7 and phase error polygons of N = ±3 are clearly visible. Our method therefore requires production of colour representations of several topography-corrected tandem interferograms with different baseline ratios (Fig. 4).
Fig. 4. Examples for (a) 12 and 13 November 1995 and (c) 31 March and 1 April 1996, of phase-unwrapping error polygons (white arrows) and (b, d) corresponding corrected interferograms.
Our proposed interactive technique to correct unwrapping errors is similar, and complementary, to the multi-baseline procedure for generating DEMs without unwrapping proposed by Reference Eineder and AdamEineder and Adam (2005). Although our manual approach is considerably more time-consuming, it is successful with small numbers of multi-baseline interferograms where an automated approach would fail.
Shifting of identified error polygons in the unwrapped topographic phase of the large-baseline differential interferogram and updating of the colour representations of the topo-corrected tandem interferograms are done in real time using ArcGIS and an Interactive Data Language (IDL)-based computer display. The values of N required to correctly raise or lower the polygons are determined by trial and error, and are iterated until no interferograms have abrupt breaks in colour (Fig. 4). Using this method, we removed all unwrapping error polygons of ∼250 pixels or larger and corrected the unwrapped topographic phase of the large-baseline differential interferogram.
To reduce the error introduced by atmospheric changes (step 4), we estimated the atmospheric contribution, φ a (Equation (1)), by using a spatial low-pass filter (cut-off wavelength ∼7 km) after first masking and interpolating all areas with motion bias or poor coherence (Reference Rabus and GhumanRabus and Ghuman, 2009). Shorter-scale atmospheric disturbances may remain; an evaluation of these remnants in the Delta River valley suggests a remaining atmospheric error of <0.2 fringes. We assume that the 1949 USGS DEM is correct in ice-free areas at scales larger than the filter wavelength. After subtracting φ a from the interferogram, we produced a topography-change map that displays the difference between the 1949 USGS data and the topographic signal in the 1995 ERS data (Fig. 5). Areas in the topography-change map that contain significant SAR layover or poor coherence were masked and populated with values from the USGS DEM. Most of these areas are in steep, north-oriented valleys outside the glacier surface. Figure 5 thus shows the actual surface change over the 46 years 1949–95, but it still contains 1949 map error. We then overlaid the topography-change map on the USGS map to produce a revised DEM, in map projection, that represents the 1995 surface of Black Rapids Glacier and its surroundings.
Fig. 5. Topography-change map, based on the difference between the InSAR DEM (1995) and the USGS DEM (1949).
Error analysis
The agreement between the ground elevations surveyed in 1995 and those extracted from the 1995 ERS DEM (see below) is evidence that the DEM we created is accurate. We quantified the error in the DEM by estimating all independent error contributions to the InSAR phase (adapted from Reference Rabus and FatlandRabus and Fatland, 2000):
where λ is the wavelength of the SAR signal (5.65 cm), r is the slant range distance (∼850 km) and θ is the incidence angle of the SAR (23.5°). Baseline errors are small for the precise ERS orbits. We estimate that, after the atmospheric correction, the combined error contribution of atmosphere, orbit and sensor-noise effects is <0.25 fringes (Reference Rabus and FatlandRabus and Fatland, 2000). The error stemming from non-constant glacier velocity is estimated from Figure 3 as <0.5 fringes. A total maximum error of ∼0.75 fringes corresponds to a possible error of ∼11 m in the ERS DEM. Dry snow is largely transparent to C-band radar, while penetration into ice is minimal (<1 m). In temperate firn, an error up to 4 m is possible (Reference Rignot, Echelmeyer and KrabillRignot and others, 2001). Our DEM, which is constructed from winter scenes with dry snow cover, therefore closely represents the ice surface in the ablation area, while in the firn area the DEM surface may be slightly biased to below the firn surface. The horizontal geolocation error is about 25 m (1.25 pixels), or 2.3 (1/tan(incidence angle)) times the elevation error.
Results
We compare the 1995 ERS DEM (Fig. 6) with several elevation datasets, including the 1949 USGS DEM, ground survey data collected from the 1970s to 1995 (M. Truffer, http://www.gi.alaska.edu/∼truffer/BRmonitoring.pdf), a high-resolution (10 m) Intermap DEM of the lower ∼10 km of the glacier in 2001, and the newly released Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEM (GDEM) based on imagery acquired between 2000 and 2009. Most of the data used to construct ASTER GDEM tiles N63W146 and N63W147, which include Black Rapids Glacier, date to before 2002 because the landslides triggered by the Denali Fault earthquake are not visible on them. We thus refer to this DEM as the 2001 ASTER GDEM. The ground survey measurements are not complete for every station. Sites that are difficult to access, such as at km 2 and km 4, were discontinued when the USGS transferred the project to the University of Alaska.
Fig. 6. 1995 ERS DEM. Glacier surface delineated by dashed boxes (a and b) is shown in more detail in Figure 7.
Comparisons with USGS DEM and ground survey data
The ERS DEM reveals much detail of the surface of Black Rapids Glacier, including lateral and medial moraines up to 35 m high emanating from the Loket tributary (Fig. 7). Pothole lakes located near the Susitna–Black Rapids divide, described by Reference Sturm and CosgroveSturm and Cosgrove (1990), are also prominent on the ERS DEM. The USGS DEM, on the other hand, shows a nearly featureless, relatively flat glacier surface with local artifacts resulting from interpolation between topographic map contours. The moraines and potholes are not evident on the DEM, although some of the moraines are visible on the original paper topographic maps.
Fig. 7. Details of DEMs (Fig. 6 for locations of (a) and (b)). Note the interpolation artifacts in the USGS DEM, the mole-run and linear artifacts in the ASTER GDEM and the moraines in the 1995 ERS DEM.
There are significant differences in the elevation of the surface of Black Rapids Glacier between 1949 and 1995 (Fig. 5), suggesting major thinning near the terminus and thickening in the accumulation area. The greatest thickening is in the upper ablation area. These results are discussed below.
Our ERS DEM matches the ground survey elevations closely, with a maximum difference of 16 m at km 29 (Fig. 8; Table 1). At all other sites where 1995 survey data exist, the ERS DEM surface differs from the survey elevations by ±3 m, which is significantly less than the calculated root-meansquare error for the ERS DEM. In comparison, the 1949 USGS DEM is typically 40–60 m lower than the 1995 survey elevations, although at the most down-glacier site (km 29) the USGS surface is 18 m above the surveyed elevation.
Table 1. Summary of changes in elevations and annual rates of change along center line of Black Rapids Glacier over different periods
Fig. 8. Comparison between the 1995 ERS DEM elevations, 1995 surveyed elevations and elevations derived from the 1949 USGS DEM, 2001 Intermap DEM and 2001 ASTER DEM. All elevations are along the center line of the glacier (see upper-right inset). Error bars represent the cumulative error of the ERS DEM and the USGS map. The inset at the bottom left shows elevations between km 5 and km 15 in more detail.
The calculated area-averaged thinning rate for the main trunk of Black Rapids Glacier and the Loket tributary between 1949 and 1995 is 0.17 m a−1. Maximum elevation changes along the center line in the ablation and accumulation zones of Black Rapids Glacier over the same period are −249 m (km 41.1) and +63 m (km 6.3). Average rates of elevation change over this 46 year period are, respectively, −5.4 and +1.4 ma−1. Ground-survey elevation changes at nearby locations are −4.9 ma−1 (km 38) and +0.5 ma−1 (km 8) for the period 1975–85, and +0.2 ma−1 (km 8) for the period 1985–94. The maximum rate of surface lowering along the center line of the glacier between 1949 and 1995 is comparable with the rate calculated from the ground survey data, which is based on the period 1975–84. The average thickening rate along the center line in the accumulation area between 1949 and 1995, on the other hand, is three times larger than the rate calculated from surveyed elevation changes between 1975 and 1985, and seven times the rate calculated between 1985 and 1994.
The glacier thickened 62 m (1.4 m a−1) in the upper ablation area (km 21), directly up-glacier of the Loket tributary between 1949 and 1995. The bulge is asymmetric along the center line: larger at the south than the north (Fig. 5). Between 1975 and 1995, the glacier has had a net elevation change of 0 m at km 20. Between 1995 and 2005, however, the glacier surface at km 20 lowered 5 m.
Comparison with Intermap Dem
The 2001 Intermap DEM extends only about 10 km up Black Rapids Glacier from the terminus, but it provides an additional check on the accuracy of the 1995 ERS DEM. Center-line elevation differences between the ERS and Intermap DEMs are small: the Intermap surface is <31 m below the ERS surface at the survey stations (Fig. 8). Ground surveys between 1995 and 2001 at km 29, on the lower part of the glacier, indicate surface lowering of ∼9 m. Both DEMs have similar surface features, and moraines are well defined on both. Down-glacier migration of the moraines over the 6 years separating the two DEMs is ∼160 ± 20 m (27 ± 3 ma−1), which is similar to the surveyed velocity of ∼30 m a−1 reported by M. Truffer (http://www.gi.alaska.edu/∼truffer/BRmonitoring.pdf) at km 29, near where the measurements were made.
Comparison with Aster Gdem
Much of the surface of Black Rapids Glacier on the 2001 ASTER GDEM is within ±10 m of the glacier surface on the 1995 ERS DEM. At the glacier terminus (km 38), the ASTER surface is 74 m lower than the ERS surface and 43 m lower than the 2001 Intermap surface. Ground surveys, however, show that the glacier surface between km 8 and km 20 changed <1 m between 1995 and 2001; in the lowermost ablation area (km 29), the glacier thinned only 10 m over this period (1.6 m a−1) (M. Truffer, http://www.gi.alaska.edu/∼truffer/BRmonitoring.pdf). The ASTER GDEM glacier surface has numerous ‘mole-run’ artifacts, which are described in the literature that accompanies the dataset (NASA Jet Propulsion Laboratory (JPL)/Ministry of Economy, Trade and Industry, Japan (METI), http://www.ersdac.or.jp/GDEM/E/image/ASTER_GDEM_Readme_Ev1.0.pdf). These artifacts are related to the linear and curvilinear boundaries that separate the various DEMs used to create the final product.
Discussion
The comparisons presented above suggest that our ERS DEM is an accurate representation of the 1995 glacier surface. The ERS–Intermap DEM elevation differences are reasonable, but apply only to the lower ablation zone. The ASTER GDEM is locally compromised by artifacts, making quantitative comparison with other datasets difficult. We therefore restrict our discussion of the glaciological implications of our work to comparisons of the ERS DEM with the 1949 USGS DEM and the ground survey data. Our discussion deals, successively, with the lower ablation zone, upper ablation zone and accumulation zone.
Lower ablation-zone thinning
The maximum thinning of 249 m in the lower ablation zone between 1949 and 1995 (5.4 m a−1), although large, is consistent with rapid wasting of a stagnant glacier terminus after a surge. Ablation of the debris-covered terminus of Black Rapids Glacier is dominated by surface erosion, rainfall and high temperatures during summer (Reference SchomackerSchomacker, 2008). Glaciokarst is prominent at the glacier terminus, a sign of rapid surface lowering. Reference Schiefer, Menounos and WheateSchiefer and others (2007) calculated thinning rates of up to 9 m a−1 between 1985 and 1999 for glaciers in the southern Coast Mountains of British Columbia, none of which exhibit surge behavior. It is difficult to extrapolate these thinning rates back to 1949 and to a different climatic zone, but they suggest hundreds of metres of thinning are possible. Reference Berthier, Schiefer, Clarke, Menounos and RémyBerthier and others (2010) report up to 6 m a−1 thinning of glaciers in the Alaska Range since 1949, which is consistent with our findings.
Upper ablation-zone thickening
Black Rapids Glacier thickened up to 62 m between 1949 and 1995 (+1.4 m a−1)over ∼6 km of its length centered near km 20 in the upper ablation zone. We offer three possible explanations for the thickening near km 20. First, it could result from an increase in mass balance due to locally higher snowfall or lower ablation. Second, it may represent dynamic thickening of the glacier due to a change in basal boundary conditions or interaction with ice from the Loket tributary. Third, it may reflect a significant error in the USGS DEM.
A local change in mass balance over the four decades is possible but unlikely. Dynamic thickening of the glacier could occur if ice flowing east from the main accumulation area of Black Rapids Glacier ponded upstream of the Loket tributary. Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996) concluded that this effect caused the glacier to thicken ∼50 m at km 20 between 1949 and 1977 (+1.8 m a−1), and surveying since then has shown further thickening of 6 m between 1980 and 1989 (+0.7 m a−1). Since 1990, however, thickening at km 20 has ceased, and M. Truffer (http://www.gi.alaska.edu/∼truffer/BRmonitoring.pdf) reported thinning of 0.2 m a−1 between 1990 and 1999 and 1.2 m a−1 between 2000 and 2009. The thickening reported by Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996) between 1980 and 1989 was preceded by 5 years of thinning of nearly the same amount, resulting in a net zero elevation change between 1975 and 1994 at km 20. Over the full ground-survey period (1975–2009) at km 20, the glacier has thinned 12 m (−0.35 m a−1). If dynamic thickening occurred at km20 between 1949 and 1975, it had to have been at a rate of ∼+2.4 m a−1, more than three times higher than the rate reported by Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996) for the 1980s (0.7 m a−1).
The behavior of the glacier over the first 35 years following the surge in 1936/37 is unknown. However, considering the typical behavior of glacier surges, it is likely that the surge thinned the glacier markedly upstream of the Loket tributary, leading to subsequent long-term thickening. The elevation of the glacier at km 20 on the 1949 USGS map is 1451 m a.s.l. The surveyed elevation of the ice surface at km 20 in 1975, 26 years later, was 1514 m a.s.l. (M. Truffer, http://www.gi.alaska.edu/∼truffer/BRmonitoring.pdf). It is important to note that the bulge reported by Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996) and shown here in Figure 5 is not the result of an increase in elevation in the 1995 ERS DEM, but rather a depression in the USGS surface (Fig. 8), which is consistent with exceptional thinning of the glacier at km20 during the 1936/37 surge.
The third explanation, that the USGS DEM is significantly in error (>60 m) in the vicinity of km 20, cannot be precluded, but it seems unlikely since photogrammetry over bare ice is generally reliable (Reference Schiefer, Menounos and WheateSchiefer and others, 2007). Figure 5 shows the shape and extent of the bulge described by Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996). The elevation changes around the bulge are smooth and gradual, further suggesting that a large photogrammetric error is unlikely. In conclusion, the bulge is most likely due to dynamic thickening of the glacier following the 1936/37 surge.
Accumulation zone thickening
The question arises as to whether the inferred thickening in the accumulation area is sensible. The inferred maximum center-line thickening at km 6.3 between 1949 and 1995 is 63 m (+1.4 m a−1). Ground surveys at the closest survey station, km 8, however, indicate only about 5 m of thickening between 1975 and 1985 (+0.5 m a−1) and 2 m of thickening between 1985 and 1994 (+0.2 m a−1). The thickening rate at km 6.3 between 1949 and 1975 would have to be nearly six times higher (+2.2 m a−1)than the rate at km 8 between 1975 and 1994 (+0.4 m a−1)for these two findings to be consistent. Given the large differences in estimated thickening rates between the repeat surveys and the 1949 and 1995 DEMs, the 1949 USGS map is likely in error in the accumulation area (Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others, 1996). Such a large error is not unexpected given the difficulties of accurately delineating the snow surface in the featureless accumulation zone of Black Rapids Glacier on aerial photographs upon which the USGS map is based.
Conclusions
We present a new method of producing accurate DEMs of alpine glaciers using ERS-1/-2 repeat-pass data and an iterative, semi-automatic, multi-baseline InSAR method. Using this method, we produce a DEM of Black Rapids Glacier that is an accurate representation of the 1995 surface of the glacier. The method can be used to create reliable DEMs of other alpine glaciers.
We use the new 1995 DEM to evaluate errors in a 1949 USGS DEM and changes in the glacier during intervening years. Maximum glacier center-line elevation changes between 1949 and 1995 are −249 m in the ablation zone and +63 m in the accumulation zone. A bulge developed in the upper ablation zone between 1949 and 1995, accompanied by center-line thickening of 62 m. Survey data collected on Black Rapids Glacier over nearly 40 years indicate that thinning over the lowest part of the glacier and thickening in the upper ablation zone are real. On the other hand, the maximum inferred thickening of 63 m likely stems from errors in the USGS DEM, although some thickening would be expected over that part of the glacier.
Acknowledgements
This research was funded through a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant to Clague, and an NSERC-PGS-D3 Scholarship, Geological Society of America Bruce ‘Biff’ Reed Research Grant, Northern Scientific Training Programme grants and an Arctic Institute of North America Grant-in-Aid to Shugar. The European Space Agency supplied ERS data (Category-1 proposal No. 3970), and M. Nolan provided the Intermap DEM. M. Truffer provided survey measurements. The ASTER GDEM is a product of the METI and NASA. Reviews by M. Braun and an anonymous reviewer strengthened the paper.
