Introduction

One of the many global effects of climate change that remain uncertain is the response of the Greenland ice sheet. The size and inaccessibility of this region, and the complexity of the system of mass discharge mean that uncertainty remains even as to the sign of Greenland’s mass balance, and thus its contribution to sea-level rise (e.g. Reference HuybrechtsHuybrechts, 1994). Some ablation estimates have been made that are based on energy-balance, degree-day and topographic modelling (e.g. Reference Braithwaite and OlesenBraithwaite and Olesen, 1990), but the more direct approach of quantifying the balance of significant outlet glaciers has not been attempted because of insufficient data.

Factors which can be used to determine the out-flux terms of mass balance for outlet glaciers include the size and characteristics of the drainage basin, the thickness of ice at the terminus, and the rate and variability of flow. Data relating to glacier thickness remain sparse, but once measured will need remeasuring only infrequently. Drainage-basin dimensions are also relatively stable and therefore can readily be derived from a variety of remote-sensing methods. Surface velocities, however, can vary on different temporal scales, depending on responses to local conditions, seasonality or climate change, and thus in situ methods are unlikely ever to yield sufficiently comprehensive data (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002).

Recent advances in synthetic aperture radar (SAR) data processing have proved effective in remotely measuring ice-surface velocities (e.g. Reference Mohr, Reeh and MadsenMohr and others, 1998). However, satellite radar interferometry (SRI) is of limited utility for measuring the velocities of major outlet glaciers in Greenland since the relevant velocities in this region are generally too high, and often only single line-of-sight measurements are available. Nevertheless, SAR imagepairs can be utilized through the recently developed technique of automatic feature tracking, based on patch intensity cross-correlation optimization, which is capable both of quantifying the relevant high glacier flow velocities ( 41m d^–1) and of simultaneously measuring both the range and azimuth components of surface displacement.

Many previous studies have applied correlation-based feature-tracking techniques to ice-surface velocity measurement using optical satellite imagery (e.g. Reference Bindschadler and ScambosBindschadler and Scambos, 1991; Reference Rignot, Kwok, Curlander and PangRignot and others, 1991; Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others, 1992; Reference DwyerDwyer, 1995). More recently, similar techniques have been applied to SAR imagery using both detected intensity correlation and phase coherence optimization (e.g. Reference Michel and RignotMichel and Rignot, 1999; Reference Gray, Short, Mattar and JezekGray and others, 2001 ;Reference WernerWerner and others, 2001; Reference Murray, Strozzi, Luckman, Pritchard and JiskootMurray and others, 2002; Reference Murray, Strozzi, Luckman, Pritchard and JiskootStrozzi and others, 2002). Data can be acquired more consistently using SAR than using optical sensors, which are subject to clouds and varying solar illumination, so monitoring ice surface velocities in this manner is much more reliable.

This paper investigates the use of SAR feature tracking to quantify the surface velocities of major outlet glaciers in a relatively large region of East Greenland.

Study Site and Data Sources

The area studied lies between the edge of the Greenland ice sheet and the east coast, between Kangerlussuaq and Scores by Sund. Glaciologically, this area comprises both glaciers draining local ice caps (e.g. the Geikie Plateau), and large outlet glaciers of the Greenland ice sheet. Of those glaciers draining locally, an unusually high proportion are known to be of surge type (Reference WeidickWeidick, 1988; Reference Jiskoot, Pedersen and MurrayJiskoot and others, 2001). This area is therefore of particular interest for studying surge behaviour and quantifying east coast ice-sheet discharge, and provides a suitable sample of glaciers for testing velocity measurement techniques on a regional basis.

The area is encompassed by 14 standard European Remote-sensing Satellite (ERS) SAR frames with slant- and ground-rangepixel size of 7.904 mby 3.962 m in range and azimuth and nominal resolution of 25–30m. Data were provided by the European Space Agency (ESA) and distributed under the VECTRA agreement by University College London. Two tandem SAR imagepairs (ERS-1and-2 separated by1day) of single look complex (SLC) images for each frame were acquired from the winter of 1995/96, making 56 images in all (Table 1). Figure 1 shows the geographic location of the region, and a mosaic of the 14 ERS-1 images. These images were chosen from the archive to satisfy the following criteria:

1. (1) Maximum coverage of the currently glacier-covered land in the region, with minimal image overlap and consecutive frames where possible.

2. (2) Minimal precipitation or periods of temperatures above zero between tandem pairs as recorded at the Scoresby Sund meteorological station.

3. (3) Minimal number of orbital frames between the tandem pairs (1cycle or 35 days for 13 out of 14 frames; 2 cycles or 70 days for the remaining frame).

4. (4) Combinations of perpendicular orbital baselines suitable for differential interferometry.

Table 1. Acquisition parameters for all images analyzed

Fig. 1. Mosaic of feature-tracking results for 35 day repeat delay at a reduced resolution (625 m×625 m pixels). Colour denotes tracked velocity magnitude (45 m d–1has maximum scale colour), while brightness indicates backscatter power for the first of the image pair. Geographic coordinates and insert show region location, while Universal Transverse Mercator (UTM) coordinates (zone 26) on 100 km square grid show scale and orientation. Coloured squares show three largest glaciers (Magga Dan in red, Christian IV in orange and Vestfjord in blue), full-resolution images of which are in Figure 2.

Using these criteria, it was possible to obtain a good coverage of the region, and Table 1 shows acquisition parameters for all images analyzed.

Methods

The underlying method of feature tracking is based on optimizing the cross-correlation between image patches in intensity-detected SAR images. Since satellite orbits are not precisely repeating, the global offset between a pair of images has first to be calculated. The difference between this global offset and the local offsets found by the peak of the cross-correlation on a regularly spaced grid of image patches then gives the displacement of the features dominating each image patch. These can be converted into a field of surface velocity estimates by taking into account the time delay between images and the slant-to-ground range geometry of the SAR imaging process. In addition, the height of the correlation peak relative to the average level of correlation may be used as a signal-to-noise-ratio (SNR) measure to filter out noisy or low-confidence measurements.

The precise method used here is described in more detail by Reference Murray, Strozzi, Luckman, Pritchard and JiskootStrozzi and others (2002) but is not significantly different from that used previously on optical images (e.g. Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others, 1992). In this study, the process of feature tracking between images, projecting into a known map reference, mosaicking the velocity measurements from many image frames and visualizing the results has been entirely automated, from the delivery of ESA SLC images to the production of geocoded mosaics of tracked features. In this way, we can be confident that these methods could be applied to other glacierized areas in a relatively straightforward manner.

A brief description of the procedure, including specific details, is as follows:

1. (1) Patch intensity correlation optimization was performed with a relatively large window size (128 by 256 SLC pixels of 7.90 m by 3.96 m in slant range) and complete overlapping coverage to find the second-order two-dimensional polynomial (of the form A+B* range +C*azimuth +D*range*azimuth, where A,B,C and D are free parameters) defining the background offset between images (the zero-velocity polynomial). Correlations with SNR 57.0 and with measured offsets more than 3 standard deviations from the fit were not used to define this polynomial, so as to minimize the effect of moving ice surfaces on this part of the analysis.

2. (2) Patch intensity correlation optimization was repeated with a smaller window size to find the patch-by-patch two-dimensional offsets between images including those due to the zero-velocity global offsets and those due to correlated surface movement. The SNRs for each patch were found simultaneously and recorded to help filter out noisy data at a later stage. Windows of 64 by 256 SLC pixels were chosen at this stage to give approximately square patches in ground range (ground-range patch size = 7.90 m 6 64/sin 23° = 1295 m; ground-azimuth patch size = 3.96m 6 256 = 1014 m).

3. (3) The background offsets were removed by subtracting the zero-velocity polynomial and displacements converted from slant range to ground range (by projection dependent on the look angle at each patch). The displacement values were scaled to m d^–1 by dividing by the delay between image pairs to create comparable tracking images for each frame of each coverage.

4. (4) The tracking and SNR images were padded to place them in a known coordinate reference with respect to a calibrated multi-looked (spatially averaged) backscatter image (10×50 multi-looking in range and azimuth, respectively).

5. (5) Finally, the two-dimensional tracking images were geocoded (using ellipsoid correction only) to a common map projection (UT M zone 26) at an appropriate pixel size (125 by 125 m) keeping track of corner coordinates as defined by precision orbital state vectors provided by ESA. Images of SNR and multi-looked backscatter were also geocoded at this stage so that they could be combined later to filter and visualize the results in a manner optimal to the whole region rather than just each single frame.

These steps were achieved for the various combinations of 1day and 35 day repeat pairs. The resulting geocoded tracking, SNR and backscatter images were then combined into geocoded mosaics of 530 by 530 km for each coverage pair.

The final stage of the process was to filter the measurements to leave only measured velocities for which there was a suitable degree of confidence. Experiments were carried out to maximize the areal coverage of measurements over glacier ice while minimizing apparent noise. Measurements were masked or rejected on the basis of four factors for which optimized values were determined:

1. (1) Application of a glacier ice mask based on analysis by Reference Jiskoot, Murray and LuckmanJiskoot and others (2003).

2. (2) Minimum SNR of 4.0.

3. (3) Minimum velocity threshold of 0.1 m d^–1.

4. (4) Maximum interquartile range of velocities within any given 5 6 5 window of 4.0 m d–¹ in either range or azimuth directions. This, along with the area-correlation nature of the process, means that there is a limit on the maximum strain rate that could be measured using this technique.

The final tracking-velocity mosaics were visualized by combining backscatter and displacement as intensity and hue to provide a combined view of surface conditions and tracked displacement. There appeared very little difference between the results from the two tandem coverages or between the two 35/70 day coverages, so only the former of each of these are discussed further.

Resolution and Error

Some previous independent estimates and measurements of surface velocities in this region have been made (Reference Olesen and ReehOlesen and Reeh, 1973; Reference DwyerDwyer 1995). However, suitably contemporary or extensive in situ measurements of ice surface velocity are not available for a full analysis of error. Measurement confidence may be estimated by analyzing the expected resolution of the technique and considering other sources of error.

With two-times oversampling, the correlation technique with the given window sizes typically gives accuracies within 1/20 of a pixel (Reference Murray, Strozzi, Luckman, Pritchard and JiskootStrozzi and others, 2002). When the range and azimuth pixel sizes are taken into account this gives a velocity resolution of 0.198 md–¹ in azimuth and 0.395 md–¹ in slant range or 1.01 md^–1 in ground range when the look angle (typically 23°) is taken into account. These values can be correspondingly reduced when features are tracked over multiple days, giving, for 35 day tracking, an azimuth resolution of 0.0057 md^–1 and a range resolution of 0.029md^–1. These orthogonal resolution estimates give a velocity-magnitude resolution of 1.03 md^–1 for 1day tracking and 0.03 md^–1 for 35 day tracking.

These values must be considered as minimum error bounds. Other sources of error include calculation of the zero-velocity polynomial which is also made through image patch correlation. Although outliers are rejected at this stage according to both SNR and global trends, the current approach does not allow moving ice (or indeed drifting sea ice) to be positively excluded from this process.

Results

Figure 1 shows the full 35day feature-tracking results. The technique has successfully mapped velocity over 57% of glacier ice in the region. It is immediately apparent that feature tracking has been most successful in measuring velocity in the trunk regions of glaciers rather than high on the plateaux and ice fields, where less coverage has been achieved. Comparison of Figure 1 with a new geographical information system-based inventory of glaciers in this region (Reference Jiskoot, Murray and LuckmanJiskoot and others, 2003) shows that the fastest measured ice flow in the region is strongly related to the drainage-basin area of each glacier. The three highlighted glaciers (Magga Dan, Christian IV and Vestfjord Gletscher) with highest measured velocities are also those with the largest total surface area (Table 2; Reference Jiskoot, Murray and LuckmanJiskoot and others, 2003).

Table 2. Matched area and maximum velocities measured the three largest glaciers in the region for 1day and 35 tracking (see Fig. 2 for location of these maxima)

Full-resolution sub-images of Magga Dan, Christian IV and Vestfjord Gletscher for the 1day feature tracking, the 35 day feature tracking and the difference between them are shown in Figure 2. The maximum measured velocity is given in Table 2, and the location of this maximum is shown in Figure 2.The coverage of results is some what better in the 1day (75–80%) than 35 day (65–75%) tracking (Table 2). However, it should be noted that in the case of 1day results the measured velocity is often lower than the resolution of the technique.

Fig. 2. The three largest glaciers in the region, showing 1day feature-tracking results, 35 day results and difference between Coverage of measured points is only reduced slightly between 1 and 35 day tracking (see Table 2), and velocities mapped different time-scales are in good agreement except at the margin. Colour scale is not cyclic (i.e. values 410 m d^–1 have scale colour). Position of maximum measured velocity is indicated by a white circle. Imaged area 25 km × 625 km.

The full-resolution images from the highlighted glaciers show the velocity typically increasing towards the centre of each glacier and within 3.5 km of the terminus (Fig. 2). However, this pattern of velocity is not always the case. For example, the northern section of the terminus of Christian IV Gletscher does not increase in velocity, whereas the southern section does. Interestingly, the maximum difference between the 35 day and 1day results also occurs close to the glacier front. Here, the velocity appears to drop very close to the glacier margin in the 1day results but not in the 35 day results.

Apart from these marginal effects, the velocities measured over 1day and 35 days are in good agreement for velocities that exceed the 1day resolution of 1 md^–1, suggesting consistency in technique.

Discussion

In general, we expect tidewater glaciers to have greater velocity than land-terminating glaciers, and that their velocity should increase towards the terminus, where buoyancy effects tend to reduce basal friction and fjord levels encourage high basal water pressure and therefore fast glacier sliding (Benn for and Reference Benn and EvansEvans, 1998).We would thus expect the regions of day strong increase in velocity close to the termini of the tidewater-terminating glaciers shown in Figure 2 to equate to regions where these effects dominate, i.e. the low-effective-pressure zone (Reference Hodgkins and DowdeswellHodgkins and Dowdeswell, 1994). However, the northerly section of the terminus of Christian IV Gletscher does not increase in velocity, and the southerly section, which does increase dramatically in velocity, moves 3–10 times faster (Fig. 2). Comparison of the velocity image with the original SAR image shows the southerly side of the glacier margin to be dominated by widespread transverse crevasses showing an extensional regime, whereas moraines on the northerly side are folded, suggesting a compressional regime. It therefore appears that only the southerly side of the glacier terminus is tidewater-terminating.

The difference image shows the marginal section of Magga Dan Gletscher to be moving significantly faster in the 1day than the 35 day measurements, while the tongue just up-glacier of this section has slowed significantly (Fig. 2). Similar features are apparent at both Christian IV and Vestfjord Gletscher (Fig. 2). Since these measurements are of ice velocity in November (1day) and through November and December (35 days), these differences are highly unlikely to be driven by surface melt or rainfall, unlike velocity increases reported at other tidewater glaciers in summer (e.g. Reference Walters and DunlapWalters and Dunlap, 1987). It is further unlikely that such a regionally consistent effect could be explained by the release of englacially or basally stored water to the terminal region. We suggest that these effects could be explained by one of two factors: (1) the annual cycle in tidewater glacier velocity, which close to the glacier margin reaches a maximum in late autumn orearly winter (e.g. Reference Van derVeenVan der Veen, 1996); or (2) the sea tide, which is known to dominate diurnal velocity variations if surface melt does not occur (Reference Walters and DunlapWalters and Dunlap, 1987), having differential effect between 1day results and those averaged over 35 days.

Conclusion

We demonstrate that intensity patch correlation optimization between 1day and 35 day repeat pairs of ERS SAR images can be used to automatically derive fields of ice-surface velocity for outlet glaciers over a large region in East Greenland.The advantage over optical techniques is that data can be collected irrespective of cloud conditions or polar night.

While 1 day tracking provides greater coverage (75– 80% areal coverage), results from 35 day tracking yield sufficient coverage (65–75%) to allow interpretation of spatial patterns of velocity over a large region of East Greenland as well as the three glaciers highlighted. The 35 day tracking results are of higher velocity resolution (0.03m d^–1) than 1day tracking (1m d^–1). In addition, 35 day repeats do not require multiple satellites to acquire, so data are much more readily available both from the archive and in the future.

The results allow interpretation of the flow patterns in large glaciers. This suggests that the SAR feature tracking may be used routinely to monitor outlet glacier velocities and temporal changes on a regional basis using any of the currently available satellite SAR sensors.