Change detection by analysis of DTMs

The terrain situation in 1984 is represented by a 10m grid DTM based on aerial photographs acquired in March 1984. Its original vertical accuracy at single pixels was limited to 35 m (Reference Schneider, Schnirch, Acuña, Casassa and KilianSchneider and others, in press b).

The second DTM represents the terrain surface situation in 2000. It is part of the freely available 90m grid Shuttle Radar Topography Mission (SRTM) dataset acquired in February 2000. Generated from the records of the C-band antenna (wavelength λ = 5.6 cm), the data subset shows a vertical error of distinctly less than 10m at each pixel (Reference RodriguezRodriguez and others, 2005). This nominal relative vertical error of 10m for the C-band is considerably less at 53˚S because vertical error in the SRTM data decreases as the latitude increases (Reference Rabus, Eineder, Roth and BamlerRabus and others, 2003).

Both DTMs were geo-referenced to a uniform 10 m grid. The subtraction of the terrain elevations of the 1984 DTM from those of the 2000 DTM resulted in a systematic spatial distribution of positive and negative elevation differences that suggested a dependency on relief exposition. As this could not be explained by natural processes, there was an obvious need to analyze this systematic error in more detail. Moreover, a two-part analysis of these elevation differences considering, on the one hand, the glaciated areas and, on the other hand, the non-glaciated areas illustrates that the ice masses of the GCN must have been subject to a considerable glacier change between 1984 and 2000. The frequency distribution of the elevation differences within the glaciated areas shows a much broader and more irregular distribution shifted to more negative values than those within the non-glaciated areas (Fig. 2). This obviously indicates that an overall loss of ice volume has occurred.

Fig. 2.

Histograms of elevation change for glaciated and non-glaciated areas of the GCN.

This finding indicates that the potential systematic error between the two DTMs is masked by additional surface changes on glacier surfaces. Only the non-glaciated areas are characterized by an almost symmetric and strongly centred frequency distribution of elevation differences (Fig. 2). This implies that the observed pattern of elevation differences between the two DTM surfaces could only have been introduced by a systematic error. The irregular changes within the glaciated areas in contrast to the non-glaciated parts must be due to real surface changes. The systematic pattern of the spatial distribution of the elevation differences between the DTMs can be very well described by computing the mean elevation differences of unchanged surfaces (e.g. rock and forest areas) as a function of exposition. The result (Fig. 3) confirms the hypothesis of a systematic error depending on exposition. This error is most probably explained by a minimal geo-referencing error resulting in a minimal displacement of the two DTMs to each other (Reference KääbKääb, 2005). Consequently, the terrain elevations of the 1984 DTM had to be corrected according to this dependency. The empirically derived deviations presented in Figure 3 were approximated by a regression in six different sections composed of both polynomial and sine fits with an explained variance of 89%. Based on this regression analysis, an exposition-dependent correction of the terrain elevations of the 1984 DTM was performed, eliminating the described exposition dependency of the elevation differences (Fig. 3). By now the result of the subtraction of the corrected 1984 DTM from the 2000 DTM showed an altitude dependency of the elevation differences. Hence, a systematic error within the SRTM data comparable with findings by Reference Berthier, Arnaud, Vincent and RémyBerthier and others (2006) must be assumed. This altitude dependency was approximated by a simple linear fit (r2 = 0:60) which was used to correct the SRTM surface. To evaluate the achieved improvement of vertical accuracy, the terrain elevations of all non-glaciated areas of the two DTMs were compared once again. The result indicates a high similarity between the two terrain surfaces. The non-corrected 1984 DTM minus the non-corrected 2000 DTM difference for non-glaciated areas shows a mean of 0.9±40 m. The difference between non-glaciated areas of the corrected DTMs amounts to –0.3±35m (Fig. 2). The given ranges indicate one standard deviation. This result is considered as a successful adaptation of the two DTMs by correcting the systematic errors inherent in both DTMs. The near-zero mean of the histogram (Fig. 2) indicates that a significant systematic elevation offset between the two DTMs after both corrections can be ruled out. Subsequently, it was possible to obtain reliable results of the volumetric glacier changes of the GCN by subtracting the corrected terrain elevations of the 1984 DTM from those of the corrected 2000 DTM.

Fig. 3.

Mean elevation changes of the non-glaciated terrain surface depending on relief exposition.

The analysis of the volumetric glacier change was performed exclusively for the 199.5 km² large GCN (Fig. 1) (Reference Schneider, Schnirch, Acuña, Casassa and KilianSchneider and others, in press b). The basic prerequisite for the computation of this volumetric glacier change is the existence of terrain elevations in both DTMs. Hence, it was necessary to omit all areas characterized by missing elevation values. Due to a lack of contrast in the underlying aerial photographs, extensive areas of the 1984 DTM surface had been interpolated. Consequently, these areas had to be treated as missing values because of the insufficient accuracy of their elevation data. In the case of the 2000 DTM, the much smaller areas of missing values are due to the specific method of terrain elevation data acquisition by means of interferometric radar remote sensing. Therefore, missing values correspond to very steep areas where shadow effects, layover and phase unwrapping problems are relevant. Thus the glaciated area where the volumetric glacier changes were derived was limited to 64.4% of the total extent of the GCN. Large gaps in data occur especially in the ice cap’s summit region (Fig. 4). To interpolate these gaps, the relation between the change in elevation of the glacier surface and absolute altitude was analyzed. By classifying the mean elevation change related to 50 m altitudinal intervals, a dependency of glacier surface change on altitude was found (Fig. 5). For its approximation, both a logarithmic and a polynomial fit were considered. The extrapolation of the logarithmic fit down to lower altitudes would have led to an obvious overestimation of ablation (Fig. 5). Therefore the polynomial fit was chosen, even though its variance of r2 = 0.90 is marginally smaller than the variance of the logarithmic fit (r2 = 0.91). From this regression analysis the estimated mean values for the change in surface elevation were obtained for the whole set of terrain elevations (Fig. 6). With the subsequent replacement of all gaps in the data, the volumetric analysis of the glacier changes of the GCN was completed. Subsequently, it was possible to derive the overall volumetric change of the ice cap and its spatial distribution. The results of this analysis serve as a basis for the constraints of possible climate forcing.

Fig. 4.

Areas of positive (DTM surface 2000 > DTM surface 1984) and negative (DTM surface 2000 < DTM surface 1984) elevation changes of the glacier surface as derived from the comparison of the two DTMs. Gaps resulting from areas of missing values in at least one of the underlying datasets are plotted as grey striped areas. Contour lines are placed at an equidistance of 250 m.

Fig. 5.

Mean elevation changes of the glacier surface depending on altitude. Data points show error bars of rms error. The plotted lines represent a logarithmic (regression 1) and a polynomial (regression 2) fit to the data.

Fig. 6.

Glacier surface change computed according to the regression analysis against altitude as presented in Figure 5. Contour lines are placed at an equidistance of 250 m.

A negative mean glacier change of 3.80m on average during the years 1984–2000 was derived for the whole GCN. This corresponds to a mean annual change of -0.35mw.e. when assuming a mean density of glacier ice of 0.9 kgm^–3 (Reference PatersonPaterson, 1994). A zone of zero elevation change at ~975ma.s.l. was identified. Below this altitude the data indicate a strongly increasing thinning of the glacier tongues, with decreasing altitude resulting in a massive loss of ice of 12.6 m on average. In contrast to that, above ~975ma.s.l. a thickening of the ice cap was observed which shows a moderate increase with altitude, resulting in a gain of ice of 7.4m on average. On an annual basis, the mean elevation change of the glacier surface approximately ranges between -2.86mw.e. at the lower ends of the outlet glaciers and +0.60mw.e. in the central part of the ice cap. Concerning the ablation areas these results are in reasonable agreement with the estimate presented by Reference Rivera, Acuña, Casassa and BownRivera and others (2002). The thickening of the central parts of the GCN, however, is unique in southernmost Patagonia to our knowledge.

Error analysis regarding DTMs

In order to check whether the elevation differences between the glacier surfaces of 1984 and 2000 are clearly distinguishable from stochastic noise, a Student’s t test was applied to the deviations between the mean differences of both DTMs in general and both the glaciated areas above and below 975 ma.s.l. It was found that the surface change of the glacier cannot be attributed to stochastic noise with a significance level of >99.9%. Consequently, the glacier change of the GCN was accepted as real change.

Error analyses shows an error of 36m for individual elevation differences between the glacier surfaces of the DTMs based on single pixels. However, the rms error over the whole domain is only 0.26 m based on the 90 m grid spacing of the original SRTM data. When each altitudinal interval of 50 m is addressed separately, the rms errors range from 1.0 m at lower altitudes to 3.9 m in the upper parts of the ice cap (Fig. 5). The polynomial fit (Fig. 5) used to approximate the altitudinal variability of the glacier change is characterized by an rms error of 5.3 m. Accordingly, the interpolated values of the glacier change show a much lower accuracy then the non-interpolated values. To assess the overall error of the mean glacier change, the rms errors of non-interpolated areas of each 50 maltitudinal interval and the rms error of the interpolated areas were combined by weighting them according to the number of pixels in each category. This calculation yields an error range of ±2.5mw.e. for the 16 year period. Hence, the annual mean glacier change of the GCN is given by -0.35±0.16mw.e. a^–1.

Reference Rignot, Echelmeyer and KrabillRignot and others (2001) report a penetration depth of the radar beam into the snowpack of 1–2 m. In the case of the GCN, where short-term melt events in summer frequently extend up to the summit area, it is reasonable to assume that this value is <1 m, so it was not taken into account in this study. Nevertheless, the penetration depth would add to the observed positive change in elevation in the accumulation area. Since the estimated precipitation offset is considered as a minimal change, this does not affect the general result obtained in this study.

Climate forcing estimated using a DDM

In order to determine the climatic forcing responsible for the observed glacier changes, a set of theoretically derived specific mass-balance gradients using varying combinations of temperature and precipitation offsets were computed by a DDM. DDMs have been studied intensely during recent decades (e.g. Reference BraithwaiteBraithwaite, 1981). Recently, reviews have been presented by Reference OhmuraOhmura (2001) and Reference HockHock (2003). The DDM calculates the ablation M at the surface, based on the positive mean daily air temperature Ti , a degree-day factor f and a stochastic term xi = 0 mmd^–1, according to

and

In this equation, f is set to 7.0mmK^–1 d^–1 for ice surface and 3.5 mmK^–1 d^–1 for snow surface according to energy-balance measurements on the ice surface of one of the outlet glaciers, Glaciar Lengua, in summer 2000 as described in detail in Reference Schneider, Kilian and GlasserSchneider and others (in press a) and on the basis of long-term measurements at ablation stakes (Table 1). The set of ablation stakes consists of 12m long, flexible, white plastic tubes that were placed into holes drilled using a steam driller. Details of the ablation stake network are presented in Reference Schneider, Kilian and GlasserSchneider and others (in press a). The measurements are in reasonable agreement with the DDM result. Differences are generally within 10%. However, at stake ‘L-up’ the deviations of 38% in 2001 and 15% in 2002 indicate a serious error which may be caused by movement of the stake within the borehole or by some other measurement error. Since the long-term measurements include periods with snow cover it is assumed that taking f for snow surfaces to 50% of f for ice surfaces is an appropriate choice, which is supported by values presented in other studies (e.g. Reference BraithwaiteBraithwaite, 1995).

Table 1.

Comparison of ablation modelled with the DDM and ablation measured at stakes

In order to consider the sensitivity of the DDM to inaccuracies of the degree-day factors, the model runs are additionally performed with f for ice surface set to 6 and 8mmK^–1 d^–1 and f for snow surface set to 3 and 4 mmK^–1 d^–1 (see Discussion).

The decrease of air temperature with altitude (–0.63 K (100 m)^–1) and the increase of precipitation with altitude (5% (100 m)^–1) were set according to measurements at automated weather stations (AWSs) at different altitudes in the years 2000–03. The details on the AWSs and the measurements for the period 2000–02 are presented in Reference Schneider, Glaser, Kilian, Santana, Butorovic and CasassaSchneider and others (2003). As a boundary condition for the start of each model run, a pattern of existing snow cover based on field observations according to Table 2 is assumed. In the course of the model run, this snow cover develops according to ablation and/or solid precipitation with respect to altitude and the meteorological record. Different snow-cover patterns at the start of each model run were considered in order to check the effect on model results (Table 2). The sensitivity of the model to this introduces an error of ±15mmw.e. a^–1. The proportion of solid precipitation, SSP, is assumed to vary depending on air temperature T according to

Table 2.

Annual mass-balance differences between DDM runs based on different starting conditions of the snow cover. Each snow pattern shows a specific variability throughout the transition zone between the ablation and accumulation zones. Snow depth is scaled linearly, using six classes between the stated upper and lower limits of altitude and snow depth. The pattern of the chosen model uses a starting snow cover of 0–500mmw.e. that is scaled between 300 and 700ma.s.l.

Equation (2) smoothly scales the proportion of solid precipitation between 100% (0˚C) and 0% (2˚C), as similarly suggested by Reference Hantel, Ehrendorfer and HaslingerHantel and others (2000) for the relation between monthly air temperature and snow-cover duration.

The DDM was extended to a surface mass-balance model for the complete GCN according to the method outlined in Reference Braithwaite and ZhangBraithwaite and Zhang (2000) on the basis of the 2000 DTM. The gaps within the 2000 DTM were replaced by terrain elevation according to the corrected 1984 DTM plus the mean glacier change between 1984 and 2000 (Fig. 6).

The DDM was run using temperature and precipitation records of an AWS named AWS GCN Puerto Bahamondes (Fig. 1) operated approximately 3.5 km from the ice margin at 28ma.s.l. since October 1999 (Reference Schneider, Glaser, Kilian, Santana, Butorovic and CasassaSchneider and others, 2003). The mean air temperature at this station between September 2000 and August 2005 is 5.8˚C, and the mean annual precipitation during the same period amounts to 4780 mm. The mass balance for this period yields a reasonable overall value of –502±15mmw.e. a^–1 according to the DDM calculation since a general retreat of the outlet glaciers can be observed. The absolute values of mass balance derived from the DDM are subject to many uncertainties and are not important in the context of this paper. However, systematic variations of precipitation and temperature yield reliable estimates of relative mass-balance changes due to climatic forcing. This makes it possible to quantify the deviation of these variables from steady-state conditions that can explain the observed surface changes. Therefore, the change in mass balance forced by temperature and precipitation offsets was computed by varying both variables stepwise. The resulting matrixes of relative mass balance of DDM runs with f over ice set to 6, 7 and 8mmK^–1 show the change in annual mass balance compared with an overall zero mass balance (Fig. 7). These data are compared with the mean annual overall elevation change of the glaciated areas derived from the DTM analysis. Both methods treat the ice cap as a whole and represent a mean glacier change including whatever thinning or thickening might have occurred within it. This means that both types of estimates represent the ice cap as a closed system including dynamic adjustment of the ice body. By interpolating the data within the matrices it was therefore possible to identify all combinations of temperature and precipitation offsets that could have led to the specific magnitudes of glacier change obtained from the DTM analysis (–346±155mmw.e. a^–1). This change is represented as bold lines in Figure 7.

Fig. 7.

Change of mass balance as a function of temperature and precipitation offsets for values of f as 6, 7 and 8mmK^–1. Changes are computed in mmw.e. a^–1 by adding constant values of temperature and precipitation changes to the input records of the DDM for the period September 2000–August 2005.