Velocity map of the Thwaites Glacier catchment, West Antarctica

Abstract

The two-dimensional surface velocity of Thwaites Glacier, West Antarctica, was mapped with 23 ascending- and 22 descending-orbit European Remote-sensing Satellite synthetic aperture radar (ERS SAR) interferograms (time range 1995–2000). The velocity map covers 175 500 km² from the Amundsen Coast to the southern turning point of the satellite orbit and comprises >80% of the Thwaites catchment. Relative velocity errors are <10% except for rare regions (about 5% of the total area) of unfavorable look geometry. Six individual tributaries were identified; their center-line velocities increase from 0 at the catchment boundary to ~0.3 km a^–1 when they join the main glacier trunk. On the main trunk, velocity increases to ~1.8 km a^–1 at the grounding line and 3.6 km a^–1 on the floating tongue. As at neighboring Pine Island Glacier, no strong longitudinal velocity gradients are found except near the grounding line. Within expected error bounds, the flow pattern appears temporally stationary, i.e. flowlines agree with the delineation of flow suggested by the pattern of velocity magnitude. A potential temporal shift of tributary boundaries must consequently be <4.4 m a^–1.

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 (Hughes, 1973; Fastook, 1984). The ice sheet’s activity is dominated by fast-flowing wet-based ice streams (MacAyeal, 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 (Hughes, 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. Bamber and Bindschadler, 1997), surface velocity (e.g. Stenoien and Bentley, 2000), grounding-line location (e.g. Rignot, 1998; Rabus and Lang, 2002) and calving events (e.g. Swithinbank 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 (Vaughan 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 (Lucchitta and Rosanova, 1997). From aerial photographs taken in 1967 (Ferrigno and others, 1993) as well as Landsat imagery from 1972 (Swithinbank, 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 (Swithinbank 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 (Lucchitta and Rosanova, 1997; Rabus and others, 2002, 2003; T. A. Scambos and others, http://nsidc.org) as well as interferometric methods (Rignot 1998). Interferometric methods have also been used to map changes of the grounding line on both glaciers (Rignot 1998, 2001; Rabus 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 (Bamber and Bindschadler, 1997) and its temporal change (Wingham and others, 1998; Shepherd 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, Stenoien 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 (Goldstein 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 (Kwok 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 (Joughin 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. Rabus 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. Kwok 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 (Bamber 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 Merson (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 (Rabus 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 (Rabus 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 (Bamber 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 (Bamber 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; Rabus 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 Vaughan 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 (Rabus 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, Stenoien 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 Stenoien 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 (Shepherd and others, 2001, 2002). Additionally, a coarse pattern of tributaries feeding the main trunk was discernible by analyzing surface topography (Bamber 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 (Stenoien 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.

5.1. Velocity field

Six tributaries feeding the main ice stream can be identified (numbered in Fig. 4). Near the margins of the catchment, ice speeds are very low (0–10 m a^–1) for most inter-tributary areas. Within the tributaries, the flow velocity increases to 35–70 m a^–1. The transition from static to moving ice is very sharp across most tributaries. An exception is tributary 3, which has no clear transition to the surrounding non-static ice (20–30 m a^–1). The flow velocity increases to about 40 m a^–1 where the tributaries join the main glacier trunk. On the tongue the velocity increases from almost 2.0 km a^–1 near the grounding line to 3.6 km a^–1 at the terminus.

As the interferometric phase on the tongue was pinned down with a feature-tracking reference measurement (October 1995–March 1996), the velocity field on the tongue represents the mean velocity during this period. Comparison of the maximum velocity on the tongue with previous investigations derived from a feature-tracking technique applied to Landsat (T. A. Scambos and others, http://nsidc.org) and ERS-1 scenes (Lucchitta and Rosanova, 1997) implies temporal acceleration. The maximum flow velocity increased by 10% since the period January 1988–January 1990 and by 6% since the period September 1993–March 1994. On the basis of ERS feature-tracking time series, Rabus and others (2003) show that the observed interannual velocity variations of Thwaites Glacier are caused by external forces (wind and ocean currents) and are restricted to the floating parts of the glacier and the area immediately inland of the grounding line.

Despite the observed variations of velocity magnitude, the qualitative spatial pattern of velocity magnitude of the Thwaites tongue has been quite stable since the beginning of observations. T. A. Scambos and others’ velocity map (http://nsidc.org) and the results of this study show similar patterns, with maximum velocities near the northwestern margin of the tongue. The interferometric derivation of the flow velocity on the tongue benefits from the invariability of the velocity pattern. By pinning the interferometric phase down to a reference value measured by the feature-tracking technique, the temporal variations of the absolute flow velocity have negligible influence on the velocity derivation.

The flow pattern on the floating tongue is remarkable: the highest velocities are observed in the western part of the tongue, while the main part of the ice mass drains through the opposite eastern section of the grounding line. This scenario is reflected by significant variations in the shear strain rates >3.6×10^–2 a^–1 (see Fig. 5c and d) and converging flowlines immediately inland of the grounding line, which results in minimum transversal strain rates <–2.5×10^–2 a^–1 (see Fig. 5b). Large crevasses formed by high longitudinal strain rates >3.75×10^–2 a^–1 (see Fig. 5a) just seaward of the western section of the grounding zone suggest that ice transport through the eastern section of the grounding line hardly compensates for the high velocities on the western tongue. High longitudinal strain rates near the grounding zone are expected as both basal friction and friction at the ice-stream margins go abruptly to zero when the main glacier trunk becomes afloat. Similar observations were made by Joughin and others (1999) for other West Antarctic ice streams.

Fig. 5 (a–c). Strain rates derived from the velocity map of Figure 4: (a) longitudinal strain rate; (b) transversal strain rate; (c) shear strain rate. Color interval is varied to provide an optimum overall impression of all three components. Some linear artifacts are caused by small steps in velocity between neighboring scenes of the velocity mosaic. (d) shows in detail the shear strain rate near the grounding zone and the corresponding subset of the velocity map. The area corresponds to the dotted rectangle in (c).

The reason for the unexpected velocity distribution across the tongue is not clear; see Rabus and others (2003) for a detailed discussion based on feature-matching time series.

5.2. Tributaries

Along- and cross-glacier velocity profiles of different tributaries are shown in Figure 6. Along-flow profiles (numbered 1–6, according to the tributaries’ labels) follow the thick flowlines in Figure 4 for each tributary. The profiles start at the inland margin of the flowlines and end at the end of the glacier tongue. The grounding-line location, derived from differential interferograms (Rabus and Lang, 2002), is indicated in profiles 1–6 as a dashed line. Two representative across-flow glacier profiles, labeled C1 and C2, are shown (denominated in Fig. 6 as profiles C1 and C2).

Fig. 6. Along-flow velocity profiles of all tributaries and representative cross-flow velocity profiles derived from the interferometric velocity map. Locations of the profiles are delineated as thick black lines in Figure 4; along-flow profiles (1–6) follow flowlines starting at the inland margins of the tributaries (km 0); vertical dashed lines indicate the grounding-line positions derived from differential interferograms. Across-flow profiles are labeled C1 and C2 beginning on the lefthand side of the tributary (looking in flow direction).

Inland of the grounding line, none of the tributaries shows a region of rapid increase in ice speed along the flow; rather, the velocity increases gradually down-glacier. The lack of rapid along-glacier velocity increase of the grounded glacier parts suggests that the ice flow is not in a bi-stable state, i.e. fast or slow. The response of the ice flow to potentially changing boundary conditions (e.g. precipitation rate, sea level, sea temperature) seems to be progressive.

A striking difference between Thwaites and Pine Island Glaciers is the degree of correlation between surface topography and tributary distribution. The main trunk of Pine Island Glacier forms within a low-slope region, surrounded by steeper slopes on three sides. All tributaries flow from higher levels into this cirque-like topographic configuration (Stenoien and Bentley, 2000). Further downstream the catchment again steepens sharply towards Pine Island Bay. The tributaries of Thwaites Glacier also join the main ice stream in a low-slope region. However, in contrast to Pine Island Glacier, the tributaries follow moderate slopes up-glacier of the confluence region (see contour levels in Fig. 3). A surface topography relationship is observed for tributaries 1, 5 and 6; the tributaries’ location is associated with a local lowering of the surface topography. A further indication for the topographic influence on the location of tributaries 5 and 6 is the existence of several nunataks and of Mount Takahe, suggesting the existence of subglacial trenches cutting through the mountain range. The corresponding velocity cross-profile is shown in Figure 6 (profile C1). For tributaries 1, 5 and 6, sliding over a hard bed is therefore more likely than sliding over a deformable bed.

The delineation of flowlines near the tributaries’ margins (see Fig. 4) is parallel to the margins of the three tributaries. Assuming a maximum error of 5° for the velocity direction (see section 4) and a flow velocity of 50 m a^–1 along the flow-lines, the temporal shift of tributary boundaries must be <4.4 m a^–1. We therefore conclude that tributaries 1, 5 and 6 are in a temporally stationary condition.

What controls the location of tributaries 2–4 in the upper central part of the drainage basin and their flow directions is less clearly defined. The velocity cross-profile through the three tributaries is shown in Figure 6 (profile C2). The transition between the tributaries is gradual, while the transition to the surrounding stagnant ice is sharp. A correlation between the surface topography (cf. contour lines in Fig. 3) and the tributary locations barely exists. Consequently, the presence of a sharp cross-flow transition in profile C2 suggests that tributaries 2–4 are thawed at the base, probably as far up as the divide. What controls the location of these tributaries and their flow directions remains unclear and is an interesting question for future investigations. Field measurements with ice-penetrating radar at Pine Island Glacier suggest that basal conditions impact the location of the tributaries there. Vaughan and others (2001) concluded that the slowly moving ice sheet is cold-based, whereas the tributaries are wet-based and lubricated. Ice-thickness data are necessary to find out whether a similar mechanism applies to Thwaites Glacier.

A potential future application may combine the along-flow velocity profiles with ice-thickness data (not acquired in sufficient detail to date) and determine the contributions of internal deformation and basal sliding (Hooke, 1998; Stenoien and Bentley, 2000). This evaluation was not done within the scope of this study, as no detailed ice-thickness data were available. The partitioning between internal deformation and sliding is an important indicator for the driving force, which controls the tributaries’ locations.

5.3. Drainage basin

The interior of the Thwaites drainage basin is not a homogeneous region passively accumulating ice and transporting ice in a spatially uniform manner to the ice stream. Rather, we find a network of fast-flowing ice tributaries interspaced by areas of nearly stagnant ice. Similar observations have been made on neighboring Pine Island Glacier by Stenoien and Bentley (2000). In contrast to the two-lobed catchment of Pine Island Glacier, the cross-section of the Thwaites catchment converges uniformly along the flow all the way to the grounding line. Vaughan and others (2001) raised the question whether the more complicated shape of the Pine Island catchment indicates unsteady conditions with respect to the generally convergent configuration typically observed for other large drainage basins. Following this argument, the convergent configuration found on Thwaites Glacier suggests a stable configuration.

Where available, the catchment margins derived by Vaughan and others (2001) are delineated in Figure 4 by dotted lines. These margins were determined by the flow-accumulation method from a DEM, starting with a fixed section of the grounding line across the ice stream and then propagating this line up-glacier in the direction of steepest ascent on the DEM. Contrary to Vaughan’s estimate, Figure 4 shows that the grounding-line section of Thwaites Glacier east must be included in the catchment. This is because both the Thwaites east ice shelf and Thwaites tongue are partially fed by the same tributary (No. 1). The flow-line pattern further suggests that a grounding-line section to the west of the tongue must be added to the catchment as well. The corresponding modified catchment margins based on the velocity field are delineated by white dashed lines in Figure 4. The area within the modified margins and the boundary of the mosaic, excluding the floating parts, is 136 000 km². By extrapolating the modified margins past the up-glacier boundary of the mosaic using Vaughan’s margins there, a drainage-basin area of 158 000 km² is obtained. As the southwestern border of the basin is probably outside the velocity map as well, this value still underestimates the actual size of the basin. Nevertheless, the measured size of the catchment already exceeds that of previous investigations, 154000 km² (Vaughan and others, 2001) and 121000 km² (McIntyre, 1984), and corresponds well to a recent measurement by Rignot and others (2002). The less pronounced correspondence between ice flow and surface topography could explain the deviating catchment sizes derived from interferometric data in this study and from DEMs in previous studies.

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.