Glacier displacement on Comfortlessbreen, Svalbard, using 2-pass differential SAR interferometry (DInSAR) with a digital elevation model

ABSTRACT Differential synthetic aperture radar interferometry (DInSAR) exploits the coherence between the phases of two or more satellite synthetic aperture radar (SAR) scenes taken from the same orbit to separate the phase contributions from topography and movement by subtracting either phase. Hence pure terrain displacement can be derived without residual height information in it, but only the component of movement in line-of-sight direction is represented in a differential interferogram. Comfortlessbreen, a recently surging glacier, flows predominantly in this direction with respect to the European Remote Sensing satellites ERS-1 and ERS-2. Four C-band SAR scenes from spring 1996 were selected because of the high coherence between the respective pairs of the 1-day repeat-pass tandem mission of the ERS sensors. 2-pass DInSAR is performed in combination with a SPOT5 (Satéllite pour l'Observation de la Terre 5) SPIRIT (SPOT5 stereoscopic survey of Polar Ice: Reference Images and Topography) digital elevation model (DEM) from 2007. The different processing steps and intermediate image products, including unwrapping and generation of displacement maps, are detailed in order to convey the DInSAR processing chain to the beginner in the field of interferometry. Maximum horizontal displacements of 18 to 20 cm d−1 in ground range direction can be detected at the glacier terminus, while a few centimetres per day characterised most of the middle and upper portions of Comfortlessbreen in spring 1996.


Introduction
Spaceborne interferometric synthetic aperture radar (In-SAR) was introduced by Goldstein and others (1988) and first applied to glaciers by Goldstein and others (1993). It takes advantage of the coherence between the phases of two SAR scenes from the same satellite orbit, even in areas with poor visual contrast like ice caps and glaciers. As it depends on active microwave sensing, SAR data can also be acquired through clouds and at night. Information on movement and elevation are both contained in a single interferometric phase (Gens and van Genderen1996;Joughin and others 1996).
Differential interferometry represents a means to separate displacement from topographic information. While Gabriel and others (1989) described DInSAR using three SAR scenes (3-pass DInSAR without DEM) over land, Kwok and Fahnestock (1996) obtained glacier velocities by differentiating the fringe structures in two independent interferograms (4-pass DInSAR without DEM). Cumming and others (1997) applied 2-pass DInSAR with DEM to measure alpine glacier flow. These different approaches have become established glacier monitoring tools, which are constantly being refined (Rosen and others, 2000;Rott, 2009). Because of the short 24 h time span and hence the high potential for coherence between SAR scene acquisitions, the ERS-1/-2 tandem mission of 1995/96 still represents a unique data source for interferometric and glaciological analyses (Eldhuset and others 2003;Wangensteen and others 2005;Rott 2009). Examples of interferometric surge research are given in Joughin and others (1996), Murray and others (2003) and Pritchard and others (2005).
This study focuses on glacier dynamics, aiming at a reconstruction of glacier flow with ERS-1/-2 data from April and May 1996 and a SPOT5 DEM from 2007. This allows a look into the past quiescent phase of Comfortlessbreen ( Fig. 1) about a decade before its surge (Sund and others 2009;Sund and Eiken 2010), to infer 1996 velocities. The goal of this paper is to render a thorough and well illustrated, low threshold description of the DInSAR workflow necessary to arrive at movement information. The entry into differential interferometry shall thus be facilitated.

Geographical setting and glacier surges
Svalbard lies in the North Atlantic (74 • to 84 • N, 10 • to 34 • E), south of the Arctic Ocean (Fig. 1). Glaciers cover about 60 % of the archipelago. Many of them are surgetype (Hagen and others 1993;Sund and others 2009) with typically low quiescent-phase velocities around 10 m a −1 (Melvold and Hagen 1998;Nuttall and others 1997). Glacier surges are sudden increases in velocity of 10 to 1000 times the quiescent flow, rapidly transferring large ice masses from higher to lower areas (Meier and Post 1969;Murray and others 2003). In Svalbard, surges tend to last in total at least ∼10 a (Dowdeswell and others 1991; Sund and others 2009) and quiescence 30 to 500 a (Dowdeswell and others 1991;Hagen and others 1993). Between surges, velocities remain too low to maintain balance, which causes build-up of glacier mass (Melvold and Hagen 1998). More knowledge on surge dynamics is important because surges can have great impact on glacier geometry, and may affect assessments relating to climate change.
Some 20 km south of Ny Ålesund on Spitsbergen, the glacier Comfortlessbreen ( Fig. 1), with an area of ∼ 65 km 2 and a length of ∼ 15 km, is flowing in a northwesterly direction from 1000 m above sea level (a.s.l.) to its partly tidewater terminus (Hagen and others 1993). Its recent surge is described by Sund and others (2009) and Sund and Eiken (2010).

Data basis
Four C-band SAR scenes from April and May 1996 from the 1-day repeat-pass tandem mission of ERS-1/-2 are used for differential SAR interferometry over Comfortlessbreen together with a SPOT5 HRS (High Resolution Stereoscopic) SPIRIT DEM from 2007 (Korona and others 2009) with 40 m resolution and 5 m height accuracy (Table 1). Each SAR image covers an area of ∼ 100 × 100 km and was provided in single look complex (SLC) format to be read into the Gamma Remote Sensing software (Wegmüller and Werner 1997) in which all data processing was done.

Differential SAR interferometry
Interferometric synthetic aperture radar (InSAR) principles An SLC SAR image s consists of amplitude, magnitude or intensity information |s| representing backscatter from the ground as perceived in Fig. 1, and phase angles ϕ recorded as fractions (0 -2π ) of a radar wavelength The phases of one SAR image alone contain few interpretable features, yet when using the phase differences or interferometric phases ϕ between two images of the same orbit and area (Fig. 2), the phase differences ϕ can be represented in an interferogram to infer ground information (Gens and van Genderen 1996;Moholdt 2010): with the distances r 1 and r 2 from a point on the ground to the satellite antenna on the first and the second pass respectively (Fig. 2) and the radar look angle θ between nadir and the target on the ground (topographic information). Those two orbit passes are separated from one another temporally by a temporal baseline B t = 24 h in the case of the ERS-1/-2 tandem mission, and in space by the spatial baseline B s , corresponding to the distance between the two satellite orbits. B s can be expressed by parallel baseline B para and perpendicular baseline B perp (Fig. 2); those are calculated and refined during interferogram generation.
In case of terrain displacement or deformation within the time interval corresponding to B t , this will lead to a distinct signal in phase difference next to the inherent topographic contribution (Fig. 2). Hence Eq. 2 needs to be expanded with a deformation term Wangensteen and others 2005): where ρ represents the displacement component in line-of-sight direction. This shows that the topography summand ϕ topo is dependent on a spatial baseline B s , while the movement summand only depends on B t , but not on any B s (Wangensteen and others 2005; Moholdt 2010).

Phase coherence and interferogram generation
Very accurate image co-registration as well as high coherence ( Fig. 3), that is only slight differences within a SLC image pair, are mandatory for interferogram generation (Rosen and others 2000;Weydahl 2001). Offset estimation via iterative image matching was used to co-register the two respective image pairs (Table 1), which were selected because of their high coherence (Weydahl 2001), with sub-pixel precision. Then one complex ERS-1 SAR image was multiplied with the complex conjugate of the other ERS-2 image taken 24 h later (Rosen and others 2000;Rott 2009). In the resulting raw interferogram, the phase differences ϕ between the two satellite images are displayed modulo 2π as colour fringes (Figs. 2, 4). Those phase differences contain a mixed signal made up not only of topography ( ϕ topo ) and displacement ( ϕ disp ) as represented in Eq. 3, but also of contributions from flat Earth trend ( ϕ f lat ), and possibly atmospheric influences ( ϕ atmo ) and noise ( ϕ noise ) (Weydahl and others 2001;Fig. 4, left): ϕ f lat and ϕ topo result from the sideward SAR viewing geometry and can be precisely computed with accurate orbit data and DEMs; ϕ disp only captures displacements in range line of sight of the sensor (Cumming and others 1989;Kwok and Fahnestock 1996;Rott 2009). Fast movements appear blurred in radar applications due to decorrelation or coherence loss (Weydahl, 2001;Eldhuset and others, 2003). Therefore, radar imagery better reflects glacial accumulation areas and slowflowing glacier parts (Goldstein et al, 1993;Joughin and others, 1996;Eldhuset and others, 2003;Pritchard and others, 2005).
Flat Earth phase removal constitutes the first step in interferometric processing (Fig. 4), followed by filtering (Rott 2009;Wegmüller and Werner 1997). Remaining phase differences are then due to topography, movement, atmosphere and noise (Fig. 4, right). When working with only one single SLC pair as in basic InSAR processing (Goldstein and others 1993), the ϕ constituents (Eq. 4) cannot be further resolved without additional supporting data; ground control points, other ground truth or DIn-SAR applications (see 4.4.) are necessary to this end. Fig. 4   adding integer multiples of 2π to ϕ whenever it jumps back to 0 from 2π ; a crucial and difficult step (Gens and van Genderen 1996;Wegmüller and Werner 1997;Moholdt 2010). Low coherence cannot be unwrapped and has to be masked out (Fig. 5). Several algorithms have been developed for unwrapping; implemented in Gamma (Wegmüller and Werner 1997) are the branchcut region growing algorithm (Goldstein and others 1988) and the minimum cost flow (MCF) technique with a triangular irregular network (TIN) (Costantini 1998). Both approaches were tested in this study.

2-pass differential InSAR (DInSAR) with DEM
DInSAR consists in subtracting the interferometric phases of one interferogram ( ϕ 1 ) from another one ( ϕ 2 ) to separate the phase contributions from  topography and movement (Eqs. 3, 4) by subtracting and thus removing either of those phases. Hence pure terrain displacement can be derived without residual height information, or alternatively pure DEM information without motion signals. The InSAR baselines B s1 and B s2 must differ for this effect, while motion is assumed to be constant (Rott 2009): A scaling factor k is introduced at this point. When aiming at the topographic phase component ϕ 1−2 = ϕ topo12 , k must equal 1 (Moholdt 2010). As topography depends on B s , while movement does not (Eq. 3), constant movement becomes fully subtracted for k = 1. Then topography or DEM information can be generated from Eq. 5. However, when aiming at displacement ϕ 1−2 = ϕ disp12 , a scaling factor k = 1 has to be used to allow for a complete removal of the perpendicular baseline and hence the topographic phase. Phase noise may increase due to interferogram scaling (Eq. 5).
Differential interferometry can be done by 2-pass DInSAR with two tandem SAR scenes and a DEM (Cumming and others 1989) as in this study, or by combining three (3-pass DInSAR) or more satellite scenes from the same orbit. Two interferograms are generated and then subtracted from one another (Joughin and others 1996;Kwok and Fahnestock 1996;Eldhuset and others 2003). In 2-pass DInSAR, a reference interferogram with ϕ corresponding to surface topography is simulated based on a DEM in radar geometry (Fig. 6). Thus the DEM must be read into radar doppler coordinates (RDC). An initial geocoding look-up table links Universal Transverse Mercator (UTM) coordinates to range-doppler geometry (Fig. 6, left). Then a SAR intensity image is simulated based on DEM and SAR imaging geometry (Fig. 6, right). Offsets are iteratively improved via pixel matching and subsequently the look-up table, so that a fine registration of the DEM becomes possible. Forward geocoding finally allows the transformation of the DEM UTM projection into RDC with improved coordinate relationships (Wegmüller and Werner 1997).
Afterwards, a first wrapped differential interferogram (Fig. 8, top) is generated by subtracting the unwrapped simulated phase k ϕ 2 (Eq. 5) from the complex, that is wrapped interferogram ϕ 1 (Fig. 4, left); both still contain ϕ f lat . Residual linear phase trends are determined by Fast Fourier Transform and removed by baseline model refinement. The corresponding unwrapped unflattened topographic interferometric phase is then resimulated with a baseline refined with help of the RDC DEM, and a second differential interferogram produced (Fig. 8, bottom).

Results
When unwrapping the interferograms, the MCF technique with a TIN (Costantini 1998) led to smoother and more continuous results than the branch-cut algorithm (Goldstein and others 1988), therefore it was given preference in this study. Fig. 7 shows the maximum spatial unwrapping extent when using 5 × 5 tiles, which divide the image into 25 rectangles. Those tiles are unwrapped individually with a certain overlap and then reassembled by the algorithm. A few horizontal and vertical linear jumps in the colour fringes remain from tiling in Fig. 7. Besides, some baseline effects from the lower right to the higher left corner can be observed (Fig. 7, left). The final interferograms were generated using 3 × 2 tiles and a reduced image subset. The elevation difference corresponding to one 2π fringe in an unwrapped interferogram is (Weydahl and others 2001):  with r a being antenna distance r 1 or r 2 . For 5-6 April 1996 with B perp ∼ 12 m, one fringe represents h 2π = 793 m of altitude, and for 10-11 May 1996 (B perp ∼ −66 m) 144 m. Glacier velocities in line-of-sight between scene acquisitions can be derived from the differential phase (Murray and others 2003;Wangensteen and others 2005) after filtering and unwrapping by MCF, leading to orthonormal displacement maps for horizontal, vertical and look vector movement component respectively (Fig. 9). The DInSAR analysis of Comfortlessbreen reveals horizontal displacements of ∼ 20 cm d −1 (5-6 April 1996) and ∼18 cm d −1 (10-11 May 1996) at the glacier terminus and < 3 cm d −1 in the middle and upper glacier portions (Fig. 9). These measured velocities correspond to <11 m a −1 along most of the glacier, while they increase towards 73 m a −1 in the lowermost 2.5 km, respectively. In the topmost area of the glacier, velocities of around 6 cm d −1 can be noted, which is commensurate to 22 m a −1 .

Discussion
The DInSAR displacement rates stated above are consistent for both April and May scenes and indicate a pre-surge velocity level on Comfortlessbreen in 1996. For Svalbard glaciers known to have surged, low velocities around 10 m a −1 are typically measured during quiescence (Nuttall and others 1997;Melvold and Hagen 1998) due to their polythermal regime. Tidewater glaciers commonly increase velocities towards the terminus (Vieli and others 2004), as is the case here (Fig. 9). The ratio from lowest (< 3 cm d −1 ) to highest (20 cm d −1 ) velocities is 7 (11 m a −1 to 73 m a −1 ). Increasing velocities in the cirque region in spring 1996 (Fig. 9, top) may represent the first indicator of surge initiation, surge stage 1 as suggested by Sund and others (2009). During surge in 2008 however, this pattern changed. Despite generally increased velocities of ∼ 2 m d −1 (730 m a −1 ), the velocity gradient, was reduced, ranging from 1.17 m d −1 (427 m a −1 ) at mid-glacier to 2.06 m d −1 (752 m a −1 ) at the terminus (Sund and Eiken 2010;M. Sund, personal communication, 2011), hence reducing the velocity ratio along glacier to < 2. DInSAR thus provides a basis for the comparison of changes from quiescence displacement levels to surge.
The potential accuracy of repeat-pass InSAR surface change detection lies in the millimeter range (Gabriel and others 1989;Goldstein and others 1993;Weydahl 2001). Yet ϕ decorrelates for ERS SAR images at ∼ 7.2 cm for horizontal movements in line of sight (Weydahl 2001). To the central west in Fig. 3, fast flowing Kronebreen with its 2 m d −1 (Weydahl 2001;Eldhuset and others 2003) appears very blue and decorrelated especially towards its margins. Such speeds as also Atmospheric influences, imprecise orbit information, signal noise and assumption of displacement along the elevation gradient represent the principal error sources in InSAR (Strozzi and others 2010). An assumed error in line-of-sight displacement of ≤ 0.7 cm for ERS-1/-2 results from a total phase error of a quarter wavelength (Crosetto and others 2008;Strozzi and others 2010). As shown in Eqs. 4 and 5, contributions from atmospheric influences ( ϕ atmo ) and noise ( ϕ noise ) are still contained in the generated (differential) interferograms (Weydahl and others 2001;Rott 2009). Ideally, ϕ atmo and ϕ noise would be fully subtracted in the DInSAR process. However, ϕ atmo may differ between repeat pass acquisitions due to varying atmospheric water vapour content. This cannot easily be corrected for, as the necessary additional data would have to cover the entirety of the atmospheric strata and is seldom available (Gens and van Genderen 1996;Hanssen 2001;Rott 2009). As shown in Fig. 3, the phase coherence of the May image set is higher than in the April one, in which some temporal decorrelation over the northern central part of the scene occurs, possibly due to snowfall, snow drift, melting or atmospheric influences (Rott 2009). Besides, SAR data penetrates the ground surface to a certain extent. Hence inferences on and even below the glacier surface are possible, concerning roughness, wetness and melting conditions amongst others.
Because of the introduction of k, ϕ noise may increase. In order to reduce noise or speckle, a salt-andpepper effect due to interferences between the many scatterers within one image pixel, the image grid can be spatially transformed to a lower resolution by averaging several looks at the same pixel, resulting in a multilook image (MLI) with better amplitude correlation and phase coherence between images at the expense of a lower spatial resolution . MLI generation represents a means to find out how much noise is left in the interferograms; therefore different MLI parameter settings should be compared to the results given here. With regard to the velocities given earlier, the errors cannot be precisely calculated within Gamma for all pixel values of a DInSAR scene. Testing all possible combinations of ERS-1/-2 SAR scenes together with more DEMs represents a means to account for the various possible error sources.
A short B perp is advantageous for deriving displacement such as glacier flow, as it leads to extremely high topographic height spans represented within one fringe, like h 2π = 793 m as given for the April pair. This reduces DEM influence to a minimum, as most fringes in the corresponding interferogram then result from movement. The tallest peaks on Svalbard reach ∼1700 m a.s.l., thus a maximum of only two entire colour cycles stem from topography in the case of B perp ∼ 12 m (Fig. 4,  top). As Comfortlessbreen reaches up to 1000 m a.s.l., some 1.3 colour cycles represent the maximum possible topographic contribution which might still be found in the unwrapped interferogram in Fig. 7 (left) which has not been differentiated from a simulated radar DEM phase. Since interferometric sensitivity to displacement is independent of the baseline B perp , while it increases with the baseline for topography, short (or zero) spatial baselines are advantageous for InSAR displacement assessments and provide better coherence (Rott 2009).

Conclusions and outlook
This study provides a step-by-step introduction to differential interferometry exemplified by Comfortlessbreen. 2-pass DInSAR with DEM yields useful and consistent displacement maps for 1996, a period where little ground truth is available for the glacier. The results allow for comparison between surge and quiescence levels, showing a reduction in velocity gradient along glacier within ∼20 years to a fourth, even though velocities increased 10 to 60 times. More flow analyses are necessary to assess surge behaviour.
It would be interesting to look further into atmospheric influences by comparing our results to local meteorological data, which are limited to ground observations in 1996. No accurate corrections are hence possible. By assessing different DInSAR combinations of the ERS-1/-2 scenes with and without DEMs, conclusions could be drawn on whether the decorrelation in the April set is due to changes in the atmosphere or on the ground.
Our results are also valuable with regard to new tandem missions which allow for interferometric analyses again. For example, the COSMO/SkyMed constellation with its current four satellites (X-band SAR) can be flown in InSAR modus with either 1-day, 4-days, 8-days or 16 days tandem recording. TerraSAR-X and TanDEM-X or the C-band Radarsat Constellation Mission represent other sources of interferometry-compliant SAR data.

Introduction
Ice shelves occupy almost 50% of the coastline of Antarctica and the majority of the continent's large outlet glaciers terminate in an ice shelf. The stability of marine ice sheets and flow rates of outlet glaciers is related to the thickness and extent of their surrounding ice shelves and margins. The collapse of ice shelves on the Antarctic Peninsula has been seen immediately to precede acceleration of tributary glaciers and thinning of ice further inland (Rott and others 2002;Scambos and others 2004). The buttressing effect of the ice shelves pinned against islands and ice rises in the Ross, Weddell and Amundsen Seas protects the unstable West Antarctic Ice Sheet from rapid collapse which could lead to a sea level rise of approximately 3.3 m (Bamber and others 2009). As well as being important for marine ice sheet stability, ice shelf velocity and thickness measurements can be used to predict inland changes by calculating variation in ice discharged across the grounding line. Variation in ice shelf extent, ice flux across the grounding line and contribution of melting and calving can be monitored by recording ice thickness and surface velocity. Four components of mass-balance are applicable to ice shelves: ice flux across the grounding line, melting/freezing at the ice-ocean interface, calving at the seaward margin and accumulation/ablation on the surface. Monitoring some of these components in combination with calculations of mass distribution, strain rate and flow patterns allows