1. Introduction

The Amundsen Sea sector of the West Antarctic ice sheet (WAIS) has been identified as being particularly vulnerable to changes in sea level and climate (Reference HughesHughes, 1973; Reference Fastook, Hansen and TakahashiFastook, 1984). The ice sheet’s activity is dominated by fast-flowing wet-based ice streams (Reference MacAyealMacAyeal, 1992). The two dominant ice streams of the region, Thwaites and Pine Island Glaciers, form the most vulnerable region of the WAIS, as they are not buttressed by substantial ice shelves. Both ice streams possess large catchments and accumulation fluxes: 154 000 km² and 52.4 Gt a^–1, respectively, for Thwaites Glacier vs 176 000 km² and 69.7 Gt a^–1 for Pine Island Glacier (D. Vaughan, unpublished information). In addition, they are largely grounded below sea level, which has been viewed as a prerequisite for unstable ice dynamics on large spatial scales (Reference HughesHughes, 1973). Consequently, repeat mapping of the spatial distribution of flow and mass balance is desirable for both ice streams. While detailed measurements of bed topography and other parameters determining basal conditions currently require a field program, other relevant parameters such as surface topography (e.g. Reference Bamber and BindschadlerBamber and Bindschadler, 1997), surface velocity (e.g. Reference Stenoien and BentleyStenoien and Bentley, 2000), grounding-line location (e.g. Reference RignotRignot, 1998; Reference Rabus and LangRabus and Lang, 2002) and calving events (e.g. Reference Swithinbank, Williams, Ferrigno, Seekins, Lucchitta and RosanovaSwithinbank and others, 1997) have already been successfully monitored with space-borne remote-sensing techniques.

Thwaites and Pine Island Glaciers enter the Amundsen Sea <200 km apart. Despite their proximity and similar size, there are striking differences between the lower reaches of the glaciers. The floating part of Pine Island Glacier is stabilized by ice shelves on both sides, which leads to more continuous calving by reducing the risk of larger calving events (Reference Vaughan, Alley and BindschadlerVaughan and others, 2001; Jet Propulsion Laboratory, California Institute of Technology, http://photojournal.jpl.nasa.gov). In contrast, Thwaites Glacier develops a large free-floating ice tongue, similar to those of Mertz and Ninnis ice streams, East Antarctica, which periodically undergo major calving events. The average velocity of the tongue, which roughly equals its growth rate during build-up, is 3.4 km a^–1 (Reference Lucchitta and RosanovaLucchitta and Rosanova, 1997). From aerial photographs taken in 1967 (Reference Ferrigno, Lucchitta, Mullins, Allison, Allen and GouldFerrigno and others, 1993) as well as Landsat imagery from 1972 (Reference SwithinbankSwithinbank, 1988) it can be concluded that the entire tongue of approximately 140 km length calved off some time in the 1940s. Remains of this iceberg, denominated B-10A or “Thwaites Iceberg Tongue”, still drift in the South Polar Sea. Between the date of the large calving event and 1988 there were a series of smaller calving events (Reference Swithinbank, Williams, Ferrigno, Seekins, Lucchitta and RosanovaSwithinbank and others, 1997), likely because the formation of a new tongue was disturbed by back-stress variations through contact with B-10A. Since 1988, when B-10A finally drifted away, the length of the tongue has increased continuously to 120 km in March 2002, when again most of the tongue calved off in a calving event (U.S. National Ice Center, http://www.natice.noaa.gov.b-22.htm) comparable to that of the 1940s.

Besides the detection of change at the calving front, data from Landsat and synthetic aperture radar (SAR) sensors have also been used to analyze the surface velocity of the lower portions of Pine Island and Thwaites Glaciers with feature-tracking (Reference Lucchitta and RosanovaLucchitta and Rosanova, 1997; Reference Rabus and LangRabus and others, 2002, 2003; T. A. Scambos and others, http://nsidc.org) as well as interferometric methods (Reference RignotRignot 1998). Interferometric methods have also been used to map changes of the grounding line on both glaciers (Reference RignotRignot 1998, 2001; Reference Rabus and LangRabus and Lang, 2002, 2003).

To determine the ice-dynamical stability of Pine Island or Thwaites Glacier, continued monitoring of the lower reaches of the ice streams must be complemented by comprehensive studies of the entire catchments. European Remote-sensing Satellite (ERS) radar altimetry has already been used to measure surface topography (Reference Bamber and BindschadlerBamber and Bindschadler, 1997) and its temporal change (Reference Wingham, Ridout, Scharroo, Arthern and ShumWingham and others, 1998; Reference Shepherd, Wingham and MansleyShepherd and others, 2001, 2002) for the Pine Island and Thwaites catchments. Over the last decade these results show a significant drop in surface elevation at the grounding lines, by about 0.75 and 2.5 m a^–1. Decrease of elevation was also found for the upper drainage basin of Thwaites Glacier, though the rate of 0.1m a^–1 is less distinct. Using SAR interferometry, Reference Stenoien and BentleyStenoien and Bentley (2000) created a large-scale velocity mosaic of the Pine Island catchment. Their mosaic reveals the spatial characteristics of the surface velocity, as well as the balance flux, of Pine Island Glacier and provides a valuable baseline for future investigations. Aiming at a similar scope for Thwaites Glacier, we have produced an interferometric velocity mosaic for the Thwaites catchment, which we present and discuss in the present paper.

2. Combination of Ascending and Descending Interferograms

An interferogram calculated from SAR repeat-pass data records surface motion within the data repeat interval (Reference Goldstein, Engelhardt, Kamb and FrolichGoldstein and others, 1993). For a non-zero interferometric baseline (which is the perpendicular component of the distance between the sensor positions of the repeat-pass data) there is an additional topographic contribution to the interferometric phase, which needs to be removed to recover the motion contribution (Reference Kwok and FahnestockKwok and Fahnestock, 1996). Smaller but unknown phase contributions are added by orbit errors or changes of atmospheric water vapor. These contributions must be included in the error budget.

Only one motion component in the satellite look direction is measured by a single interferogram. A possible way to estimate the motion of the Earth’s surface in three dimensions is to use the motion phase contributions from interferograms of ascending and descending look geometry plus an additional assumption on the motion direction (Reference Joughin, Kwok and FahnestockJoughin and others, 1998). According to Figure 1, the solution for the velocity vector v = (υ_x, υ_y, υ_z ) is implicitly given by the equation system

Fig. 1. Velocity derivation from ascending and descending interferograms. v: surface velocity vector; n: normal vector (0,0,1) on horizontal plane; v · n _asc/desc: velocity vector components in satellite look directions; θ_asc/desc: look angles; γ_asc/desc: satellite heading angles with respect to the x axis; γ₀: angle of flow direction.

Here, n is the unit vector in ascending (asc) and descending (desc) look direction, respectively. υ _look = Φ_look · (4π · Δt) denotes the measured velocity components in the satellite look direction, with Φ_look being the corresponding interferometric motion phases and Δt being the time interval between the first and second image acquisition. The solution to the equation system (1) is given by

One possible choice for n is in the xy-plane and perpendicular to the horizontal flow direction. Horizontal flow direction is locally approximated by moraine features or the down-glacier gradient of a digital elevation model (DEM). This choice allows for the derivation of the vertical velocity component (e.g. Reference Rabus and LangRabus and Lang, 2000). Alternatively n can be chosen as the surface-normal vector, which is equivalent to assuming surface-parallel flow. For the derivation of the velocity field of Thwaites Glacier, the ice surface is approximated as a horizontal plane. In this case n = (0,0,1) is constant for each pixel and we obtain

for the solution to Equation (1). Here θ and γ respectively denote satellite look and heading angles as defined in Figure 1. The application of Equation (3) for regions of the mosaic that have ascending and descending coverage yields magnitude and direction of the surface velocity.

In regions with single interferometric coverage, flow features visible in the image amplitudes and the inter-ferogram are digitized to determine the horizontal direction of the ice flow. The underlying assumption of stationary flow, together with digitization and interpolation errors, leads to a higher error budget for the single coverage regions. The projection of the measured velocity components onto the given flow direction is defined by (symbols refer to Figure 1):

In order to produce a homogeneous flow field in the single coverage regions and to minimize triangulation errors between digitized vectors, flow features were chosen area-wide at coarse spatial separation of a few kilometers.

3. Data Processing

The extent of the partial mosaics from ERS frames of ascending and descending orbits is shown in Figure 2a and b, respectively. Several processing steps must be carried out before applying Equation (3) to the unwrapped ascending and descending phase mosaics. The topographic-phase contribution in each interferogram was eliminated by subtracting a corresponding topography-only interferogram. The best way to obtain a topography-only interferogram is by subtracting two interferograms of the same look geometry but different baselines (e.g. Reference Kwok and FahnestockKwok and Fahnestock, 1996). However, as no multiple interferometric coverage of the same look geometry was available at the archive of the German receiving station O’Higgins, we instead calculated a synthetic topography-only interferogram from a coarse Antarctic altimeter DEM (Reference Bamber and BindschadlerBamber and Bindschadler, 1997). The vertical accuracy of this DEM is ~1.4 m in flat regions, and ~10.7 m in steeper regions with slopes up to 1.1%. The mean slope of the mosaic area is 0.3% which consequently allows an average topographic elimination within 0.12 fringes in the center of the scene, with the maximum baseline used for the mosaic B _eff = 293 m. For smaller effective baselines, the accuracy of topographic elimination is proportionally increased. Due to the 5 km pixel spacing, small-scale surface undulations are not represented by the topographic correction. Hence, the derived velocity field appears locally non-uniform in some slow-moving areas.

Fig. 2. Locations and effective baselines of ERS frames plotted on the Advanced Very High Resolution Radiometer (AVHRR) Antarctica mosaic of Reference MersonMerson (1989) (projection is polar stereographic). (a) Ascending mosaic; (b) descending mosaic.

The mosaicking of unwrapped, topographically corrected and geocoded interferograms requires consistent calculation of the unknown phase constant for each interferogram. Common pixels of N overlapping scenes are linked via the relation

that describes the transformation from slant range to horizontal motion phase. Here, denotes the motion phase projected onto the horizontal plain, and i (1 ≤ i ≤ N) is the overlapping interferogram number. The const _i are the unknown phase biases, are the corresponding motion-phase values, and θ ⁽ⁱ⁾ are the corresponding incidence angles of pixels in the overlap area (incidence and look angle are identical for the horizontal-flow assumption). For two neighboring scenes of the same orbit, the condition θ ⁽¹⁾ = θ ⁽²⁾ requires that one of the phase biases be given; the calculation of the difference const₂ — const₁ is then trivial for this case. For the more general case of overlapping scenes from different orbits (of same look direction, either ascending or descending), the dependency on θ should theoretically allow us to find both biases simultaneously. A least-squares fit of Equation (5) over all pixels of an overlap area that resulted in solutions for both const₁ and const₂, however, turned out to be unreliable. This is apparently due to baseline errors causing uncompensated phase ramps in some of the scenes (see corresponding discussion in section 4). Therefore, in a simpler approach, one fit parameter const₁ or const₂ was eliminated by defining a reference scene with known phase bias. A feature with known velocity (e.g. a non-moving nunatak) was identified in one ascending as well as one descending scene (location marked as a cross in Fig. 4). The phase values at the corresponding pixel locations were then subtracted as biases from the reference interferograms. Both interferograms are then transformed to horizontal motion phase with Equation (5) and remain unchanged throughout the following mosaicking process.

Starting with interferograms overlapping those of the reference orbits, all interferograms with similar look directions were matched successively. Equation (5) was simply solved by fitting the remaining bias plus an additional linear two-dimensional-phase ramp in the least-squares sense within each overlap area. Phase ramps induced by baseline errors are eliminated to first order by this method. The phase values were projected to the horizontal plain by subtracting the fitted ramp and bias. The procedure was repeated until all overlapping interferograms had been assigned a phase bias and a phase ramp. In the descending mosaic a gap with no data separates the two southernmost orbits from the others. The gap’s width ranges between 20 and 40 km (see Fig. 2b). Here, automatic fitting was not possible. Ascending interferograms and neighboring descending interferograms indicate that the interferometric phase varies only slightly in that area. Hence, values for bias and ramp were empirically determined in order to equalize the phase values at opposing margins of the two mosaic parts. Finally, the surface velocity was evaluated with Equation (3). Temporally variable tidal uplift causes a phase bias on floating regions and a vertical motion component. Both effects significantly bias the interferometric phase seawards of the grounding line. The Thwaites Glacier tongue also shows a time-variable rotation (Reference Rabus and LangRabus and others, 2003), which additionally impacts the phase and has to be taken into account when calculating the horizontal surface velocity. A simple approach was used to derive the correct motion field on three different floating regions: the glacier tongue and the adjacent regions east and west of the tongue. The regions were delineated by the grounding line, the location of which was determined from differential interfero-grams (Reference Rabus and LangRabus and Lang, 2002). Using the feature-tracking technique, velocities at reliable tie points were determined for each region. The reversal of Equation (3) delivers the unknown biases, which were added to the measured ascending and descending components υ _asc and υ _desc. Within the individual regions the biases between interferometry and feature-tracking velocities were assumed to be constant.

4. Error Analysis

An error budget was estimated for the magnitude and direction of the velocity vector. There are four potential contributions to the error budget: (i) non-horizontal flow, which results in an erroneous choice of n in Equation (1); (ii) incomplete topographic-phase removal; (iii) incomplete removal of phase biases on floating glacier parts; and (iv) errors associated with single interferometric coverage.

1. (i) The error associated with the simplification of horizontal flow is estimated from Equation (2) If the additional plane containing the velocity vector is sloped (locally) with respect to the assumed horizontal plain represented by n = (0, 0, 1), the normal vector of the true plain can be written as n + Δn. The corresponding residual velocity error vector Δv can then be derived from Equation (2) to first order as

   

   In coastal regions, ascending and descending orbits cross at angles close to 90°, which minimizes numerical errors associated with a small denominator n · (n _asc × n _desc). Near the turning point of the orbit at the southern margin of the mosaic, the relative ascending- and descending-orbit orientation is unfavorable, as both flight directions are almost parallel. Here, the influence is more significant. The slope of the glacier surface and the related error Δn were derived from an Antarctica DEM (Reference Bamber and BindschadlerBamber and Bindschadler, 1997). Due to the coarse resolution of this DEM (5 km), the error map represents a spatial average. A relative-error map of velocity magnitude Δυ = |Δv| is shown in Figure 3 for the entire area of the mosaic. A contour map of the ice surface topography is superimposed.

Fig. 3. Map of the relative error Δυ/υ of the interferometric velocity field referred to the horizontal-flow assumption. Superimposed are contour levels of the surface topography derived from an altimeter DEM (Reference Bamber and BindschadlerBamber and Bindschadler, 1997), smoothed by 40 km × 40 km filter. Contour interval is 100 m.

1. For the floating glacier parts with virtually no relief, the error contribution from Equation (6) is expected to be negligible. Topographic relief is also generally low for most of the area inland of the grounding line, reducing the relative error to <2% there. This even holds for the southernmost parts of the mosaic, where flight directions of ascending and descending orbits converge. An error of 2–5% is expected for coastal areas above the grounding line. These areas comprise the steeply inclined margins of the Antarctic ice sheet (see contour lines in Fig. 3). To reduce the contribution of the topography to the velocity error, a surface-parallel-flow instead of a horizontal-flow assumption could be used. However, the coarse resolution of the DEM, which we need to determine the surface slope, as well as the DEM’s unreliability in steeper regions, would cause other errors instead. The creation of a higher-resolution interferometric DEM is not possible due to lack of multiple scene coverage in the archive.

2. (ii) The second error contribution due to the incomplete topographic-phase removal is included as a 0.2 fringe background in the error map weighted by the effective baseline. This value was chosen to account for the steepest areas of the mosaic and local topography undulation not resolved by the DEM.

3. (iii) Errors on the floating glacier parts due to incomplete removal of the phase biases there are associated with the accuracy of the feature-tracking technique (about 50 m a^–1; Reference Rabus and LangRabus and Lang, 2003). As the surface is almost entirely flat and the flow direction is similar on all floating parts, the error is assumed to be constant there. The error map contains a constant contribution of 50 m a^–1 on the floating parts, neglecting the higher error within the spatially limited width of the flexing zone along the grounding line. Here, errors of several fringes are possible, diminishing within a few kilometers away from the grounding line to the estimated 50 m a^–1.

4. (iv) An accuracy of 5° is assumed for the flow direction derived from manual digitization of flow features and subsequent triangulation. For areas with single interferometric coverage this translates via Equation (4) into a corresponding velocity error that is included in the error map of Figure 3.

   In summary, the error of the velocity magnitude is small (<10%) for about 95% of the total area of the mosaic. In slow-moving areas error contributions (i) and (ii) (to a minor degree) result in relative errors >10% though absolute errors are very small here. The large relative error for areas of single coverage is the estimated maximum error. The velocity direction derived for the ascending/ descending coverage guided the digitization of the flow features in the adjoining areas of single coverage. Consequently, the error is smaller near the boundary of the ascending/descending coverage, and the velocity map (Fig. 4) also shows no step across the boundary.

Fig. 4. Interferometrically derived velocity map of the Thwaites Glacier catchment superimposed on an amplitude mosaic. The colors indicate the velocity magnitude; flow directions are delineated by black flowlines. Thick flowlines delineate the location of the profiles shown in Figure 5. The white cross indicates the location assumed to have zero velocity. Numbers are tributary identifiers. Dotted and dashed white lines delineate the margins of the drainage basin derived by Reference Vaughan, Alley and BindschadlerVaughan and others (2001) from an altimeter DEM and from this velocity map. The white circle on the floating tongue marks the location of a bull’s-eye pattern in the interferograms indicating a permanently grounded zone (Reference Rabus and LangRabus and others, 2003). The grounding-line location is indicated as a black line.

Of further interest is the accuracy of the derived motion direction represented by the heading of flowlines shown in Figure 4. The error of the velocity direction Δα_v is

As with the error map of the velocity magnitude shown in Figure 3, we used Equation (7) to calculate an error map of the velocity direction. Evaluating Equation (7) with Δv taken from Equation (6), we find that the contribution to the angle error between the calculated and actual velocity direction caused by the horizontal-flow assumption is negligible for most parts of the mosaic, and it is therefore not shown. Only the projection of Δα_v into the xy-plane is considered. The error contribution is generally <0.5° and reaches a maximum of 2° for the steepest regions. The error is reduced in the steeper coastal regions by ascending and descending orbits having similar orientations towards the flow direction there. Error contributions (ii) and (iii) are negligible, as potential flow-angle errors will be significantly less than 1°. The largest impact on the accuracy of the ice-flow direction is expected from contribution (iv). Possible errors come from inaccurate digitization of the flow direction in areas without clear flow features and from subsequent triangulation of the digitized flow vectors. A maximum error of 5° is assumed for contribution (iv).

Within their study of the neighboring Pine Island catchment with ERS SAR interferometry, Reference Stenoien and BentleyStenoien and Bentley (2000) had to assume baseline errors as large as 80 cm. They concluded that the precise baseline product for ERS is significantly less accurate for areas close to the southern turning point of the orbit (compared to nominal errors of ~8 cm). In the present study, phase ramps between overlapping scenes of neighboring orbits were observed mainly in the range direction for several frames near the southern orbit turning point. These phase ramps could be caused by incomplete elimination of the topographic phase (included in Fig. 3) and/or incomplete elimination of the “flat-Earth” contribution due to baseline inaccuracies. The latter error contribution is not included separately in Figure 3.

Following Reference Stenoien and BentleyStenoien and Bentley (2000), the overall accuracy of the derived surface velocity is estimated by comparing the measured flow velocity at presumably stationary points in opposite parts of the mosaic. Flow velocities at these points, which lie between two ice-stream tributaries or near rock outcrops, are on the order of 0–10 m a^–1. Therefore, for slow-moving regions an accuracy of 10 m a^–1 is assumed. Locally, small-scale surface undulations which are not resolved by the elimination of the topographic phase cause larger errors.

While the ascending scenes were all acquired within 4 months (November 1995–March 1996), descending scenes from different years (February 1996–February 2000) were used for the mosaic under the premise of a maximum spatial coverage. In another study (Rabus and others, 2003) it was found that velocity variations >30 km above the grounding line were very small (<25 m a^–1) during the observation period. No seasonal variations were found. In contrast, the mean and maximum velocity on the floating glacier tongue underwent significant variations of about 10% of the total velocity during the observation period. Potential influence of the variations on the velocity map is minimized by pinning both the ascending and descending phase mosaic down to a feature-tracking reference value derived within the period of the ascending acquisitions. The error resulting from temporal velocity variations is therefore not included in the error map.

5. Velocity Map and Interpretation

Previous studies of Thwaites Glacier delineated the downstream 33 000 km² of the main glacier trunk (Reference Shepherd, Wingham and MansleyShepherd and others, 2001, 2002). Additionally, a coarse pattern of tributaries feeding the main trunk was discernible by analyzing surface topography (Reference Bamber, Vaughan and JoughinBamber and others, 2000). Details of the velocity field of the upper Thwaites catchment and its tributaries, however, were unknown. In this paper, we mapped >80% of the glacier’s catchment. The precise delineation and detailed velocity distribution within the network of tributaries that were assumed to feed the main trunk of Thwaites Glacier are presented. Similar tributaries have been found to feed Pine Island Glacier (Reference Stenoien and BentleyStenoien and Bentley, 2000).

Figure 4 shows the interferometrically derived velocity map of the Thwaites catchment. Velocity magnitude is shown with a color bar, and flow directions are delineated by flowlines derived from the velocity field. For better orientation, a mosaic compiled from ascending and descending ERS backscatter amplitudes is underlain. Rock outcrops and areas of single interferometric coverage where ice flow is almost perpendicular to the satellite look direction were excluded from the velocity derivation and are masked out in the map. The reference point assumed to have zero velocity is located on the foot of the southern edge of the volcano Mount Takahe (3460 m); it is marked with a white cross in Figure 4.

6. Conclusion

The presented velocity map provides an accurate benchmark for future investigations of Thwaites Glacier and its drainage basin. We have identified six major tributaries of Thwaites Glacier and revealed several interesting details of the surface velocity field and longitudinal, transversal and shear strain patterns. The flow and location of three tributaries is likely to be controlled by subglacial topography. No indication could be found that the latter is also the dominant mechanism controlling the location of the other tributaries. The fact that surface flowlines appear largely parallel to the ice-stream margins (the latter are clearly visible through the sharp contrast in velocity magnitude there) indicates to us that the configuration of the tributaries may be stable on long time-scales. However, many of the other observed surface velocity features, such as the conspicuous bending of ice flow-lines in the main trunk about 50 km upstream of the grounding line, cannot be safely interpreted without detailed ice-thickness measurements, which are not available at present. On the other hand, knowledge of the precise location of these surface velocity features will greatly assist the planning of future field transects with ice-penetrating radar or other geophysical methods.