Balance velocities (U _b)

Here we have used a series of recently published datasets, derived from satellite remote sensing and in situ measurements, to calculate a quantity known as the balance velocity (U _b). This quantity is derived by calculating a balance flux and dividing this by the ice thickness. It was calculated using the approach proposed by Reference Budd and WarnerBudd and Warner (1996), which is outlined below.

If A is the net surface mass balance in ice equivalent units per unit area per unit time then at any point with ice thickness Z and depth-averaged horizontal velocity U _b, steady state is obtained if:

U _b is consequently the horizontal vector flux density in the direction of flow divided by the ice thickness:

where Φ(x, y) is the balance flux distribution. For any closed area, S (in this case the drainage basins of PIG and TWG), the total gain of mass through the upper or lower surface is balanced by the net outward flow through the boundary:

where ∂l is an elemental length along the boundary C of the surface, n is the normal to the boundary and Ψ is the scalar flux integral (or volume flux) across the boundary C. Here a finite-difference scheme was used on a quadrilateral grid allowing the balance flux distribution to be separated into its x and y components based on the slope direction, θ:

where

The computer scheme used to calculate the scalar fluxes has been presented in detail elsewhere (Reference Budd and WarnerBudd and Warner, 1996), and only the most relevant points are considered here. The key to the scheme is that rather than calculating vector flux distribution, Φ, the total scalar flux, Ψ, is calculated for each gridpoint by discretizing Equation (3) in the form:

where Ψ^(out) and Ψ⁽ⁱⁿ⁾ are the outgoing and incoming total scalar fluxes respectively.

A centred difference approach was used to estimate slope and flux vectors, which means that for a point x _i,j only the four surrounding points, x_{i −1,j} x_{i ,j−i} x _i+l,j x _{i,j +1} are relevant to the properties of x _i,j . Thus flux from x _i,j can only flow into or out of one of these four points, and Ψ⁽ⁱⁿ⁾ is determined by the slope vector calculated using the four surrounding points and the magnitude of Ψ^(out) in the upstream gridpoints. Thus, in principle, flux should be conserved. This is the case except in the vicinity of bedrock outcrops, where the scheme does not always avoid flow into an area of zero ice thickness. This is, however, a minor effect, and comparison of the net mass balance, A, with the flux crossing the boundary of S indicates that flux is conserved to better than 0.1%. The key to the scheme is the sorting of gridpoints in order of descending elevation, so that the total flux into a point is known from the contribution of upslope points that have already been processed. This allows for a simple calculation of the flux integral leaving any gridpoint x_{i ,j} . For points with equal elevation, their sorted order depends on where they lie, spatially, in the two-dimensional array. To ensure conservation of flux, it is important that there should be no elevation hollows, where flux can flow in but not out. To avoid this, elevation “polishing” is carried out. This involves identifying and smoothing hollows by adding a small elevation to the points within the hollow. This procedure does not affect the pattern of flow but ensures conservation of flux (Reference Budd and WarnerBudd and Warner, 1996).

The accuracy of the calculated balance velocity is dependent upon the quality of the input data and the method used for calculating them. Equations (1) and (2) indicate that errors in U _b are proportional to errors in accumulation and inversely proportional to errors in ice thickness. Errors in slope aspect have a regional impact as they affect the area, S, over which the integration of A takes place, and there is therefore no direct relationship between these errors and the uncertainty in U _b. A reasonable approximation, however, is that errors in area scale linearly with errors in U _b. Slope aspect errors also have a local effect on the convergence or divergence of flux at any gridpoint, and can therefore influence the local pattern of flow even if the flux, integrated over the entire area, is correct.

Balance velocities reproduce the present-day, depth-averaged velocity if, and only if, (i) there are no errors in the input datasets, (ii) the drainage basin is in balance, (iii) there is time-invariant flow and (iv) the ice flow obeys the dynamical assumptions used in the calculation of U _b. Here, we will consider the role of each of these factors in determining the reliability of the balance velocities used in this study.

Input datasets

Three input datasets are required to estimate U _b: surface slope, ice thickness and mean net surface mass balance. A DEM, containing information on ice surface slope, was derived from ERS-1 satellite radar altimetry and terrestrial data (Reference Bamber and BindschadlerBamber and Bindschadler, 1997). For slopes less than about 1°, the bias in the DEM is typically <10 m (Reference Bamber, Ekholm and KrabillBamber and others, 1998). This equates to an error in slope, over 5 km, of about 10%. The slopes for TWG are <1° (Fig. 2a). This is also the case for the central part of PIG, although there are higher slope areas (up to 1.5°) near the boundaries of the glacier (Fig. 2a). The accuracy of the DEM in the area of interest was examined using surface elevation data obtained from two of the RES flight-lines shown in Figure 1. The positions of the flight-lines have been plotted on a shaded relief plot of surface elevation (Fig. 2b) to indicate their position in relation to the local topography. The data for PIG were collected in 1998 using kinematic global positioning system for navigation control and should be accurate to ±5 m. The data for TWG were flown during the 1970s, and the accuracy of the navigation for these data is considerably poorer. Accuracy in elevation decreases with distance inland, and is typically on the order of several tens of metres. Figure 3a shows the comparison between the DEM and the airborne elevation measurements on PIG. There is good agreement for the floating tongue, but the break in slope in the DEM precedes the position indicated by the airborne data by about 10 km. This flight-line lies close to the lefthand margin of PIG where the surface slope is at its steepest (Fig. 2b). Figure 3b shows the comparison for TWG, and it is apparent that there is relatively good agreement along the length of the profile, which crosses the glacier perpendicular to the flow direction as well as down-glacier (Fig. 2b).

Fig 2. (a) Surface slopes of 0–1° derived from a 5 km spacing DEM of the ice sheet (Reference Bamber and BindschadlerBamber and Bindschadler, 1997). (b) Shaded relief plot of surface elevation covering PIG and TWG. The locations of the elevation profiles plotted in Figure 3 are indicated by the thick solid lines.

Fig. 3. (a) Comparison of elevations on PIG. The solid line is from airborne RES measurements, and the dashed line shows elevations interpolated from the DEM. (b) As for (a) but for TWG. The locations of both profiles are shown in Figures 2b, 4 and 5.

Errors in elevation are spatially variable and are difficult to relate directly to errors in U _b. However, based on the lower slopes around TWG (Fig. 2a) and the comparison shown in Figure 3, we believe the accuracy of the DEM, in the vicinity of the fast-flowing ice, is better for TWG than for PIG. It is probably comparable for the inland catchment areas, where the slopes are similar and low in both basins. We have assumed that the uncertainty in the area, A, over which the integration in Equation (3) takes place is 10%. This is based on a comparison between the area of the PIG drainage basin derived from the DEM used here and a different one (Reference Vaughan, Alley and BindschadlerVaughan and others, 2001).

The presence of longitudinal stresses in the ice can result in local changes to the ice-flow direction. Theory suggests, however, that over distances of 10–20 Z, the effect of longitudinal stresses is averaged out, and it can be assumed that ice flows in the direction of the steepest slope (Reference PatersonPaterson, 1994). This assumption is made in the calculation of U _b ¹, as discussed earlier. As a consequence, it was necessary to smooth the DEM to reduce the influence of longitudinal stresses on the flow direction. The effect of the smoothing distance on the balance velocities is to redistribute flow on scales similar to the smoothing distance (see, e.g., Reference Bamber, Hardy and JoughinBamber and others, 2000a, fig. 8). The DEM was smoothed over a spatially variable distance of 20 Z. This value was found to provide the best agreement between U _b and InSAR-derived velocities in Greenland (Reference Bamber, Hardy and JoughinBamber and others, 2000a). The true scale over which longitudinal stresses can be ignored, however, is a function of the ice rheology and basal traction, which are spatially variable quantities. There is, therefore, no single length scale or factor of ice thickness that will be correct everywhere.

The ice-thickness grid was obtained from a new compilation of data known as BEDMAP (Reference Lythe and VaughanLythe and others, 2001) and is an improvement on the data used to calculate balance velocities previously published for Antarctica (Reference Budd and WarnerBudd and Warner, 1996; Reference Bamber, Vaughan and JoughinBamber and others, 2000b; Reference Huybrechts, Steinhage, Wilhelms and BamberHuybrechts and others, 2000). Ice-thickness measurements were obtained from airborne RES surveys, supplemented with ground-based spot measurements from seismic soundings (Fig. 1). Both glaciers have RES flight-lines along their length and at least one cross-profile. Ice thickness is also relatively well constrained near the grounding line, based on an inversion of surface elevation, assuming hydrostatic equilibrium (Reference Lythe and VaughanLythe and others, 2001). The accuracy of Z calculated in this way is about 20%. The RES data have an accuracy within the range 1–10%, dependent primarily on the accuracy of the navigation. This, in turn, is dependent on when the data were collected, with data from the 1970s having the poorest navigational accuracy. For PIG, the RES coverage comprises a mixture of high-accuracy, recently acquired data and older, less accurate data. For TWG, coverage is provided by 1970s data only (Reference Lythe and VaughanLythe and others, 2001).

It is not only the quality of the raw data, but also its spatial density, that determines the accuracy of the gridded dataset. The accuracy of gridpoints is therefore dependent on their proximity to measurements. Points that are distant from any measurements could have relatively large errors. For the areas of relevance here (in and around the faster-flowing parts of the glaciers) the maximum distance to a measured or estimated ice thickness is about 30–40 km. The mean bedrock slope and ice thickness for the area covered by the interferometry are 1.4° and 1205 m, respectively. A trough or peak with this mean slope and a width of 30 km would have an amplitude of 367 m, which is 30% of the mean thickness. The rms of this and the random error in each thickness measurement is 32%. The error in Z varies, therefore, between about 10% and 32% depending on the proximity to a measurement. We have used the average of these two limits (21%) as an approximation for the mean error in Z over the area covered by the interferometry.

The grid of mean net surface mass balance was obtained from a recent compilation that combined in situ and passive-microwave satellite measurements (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999). Uncertainty in the surface mass balance over the whole ice sheet is estimated to be 5% (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999). Much of the ice-core data used to estimate accumulation rates lies in West Antarctica, with a relatively high concentration of points in the Pine Island and Thwaites basins (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others 1999, fig. 4). There is, however, a strong gradient in accumulation rate between the coast and inland areas, particularly for the Pine Island basin. For this basin, therefore, we have assumed that the uncertainty in A is 10%. For Thwaites, there are fewer in situ data (Fig. 1) but still substantially more than the average for the continent as a whole. As a consequence, we have assumed that the uncertainty in A for the Thwaites basin is 15%.

Fig. 4. Plots of (a) InSAR velocity (U _insar ) and (b) balance velocity (U _b ) for PIG. The direction and magnitude of flow is indicated by the black arrows, with the maximum length equating to a velocity of 1000 m a⁻¹ The grounding line, as determined by interferometry and used to estimate fluxes, is shown as a solid blue line. The locations of ice-thickness flight-lines are shown by the green lines. The dashed sections indicate the locations of the elevation profiles plotted in Figure 3. Ice thickness was also determined close to the grounding line, based on the assumption of hydrostatic equilibrium and a knowledge of the surface elevation and densities of sea water and ice. The black areas within the area of coverage of the SAR data are where it was not possible to obtain interferometrically derived velocity.

Based on the uncertainties in the three input datasets discussed above, the rms error in U _b for PIG and TWG, due solely to errors in the input data, is 25% and 28%, respectively. It should be emphasized that these are mean values for the areas of interest and that there is considerable spatial variability about these means.

Steady-state flow assumption and time constants

U _b is a time-averaged quantity. The approximate time constant for the response of a drainage basin in Antarctica is typically on the order of 10 000–20 000 years (Reference PatersonPaterson, 1994). This is approximately the time it would take the topography to relax to a new equilibrium condition given an instantaneous change in boundary conditions or inputs. Locally, the relaxation time is inversely proportional to the velocity. Thus, faster-moving flow features within a basin (such as an ice stream) will have a considerably shorter time constant than the average for an entire basin. Balance velocities are therefore representative of the depth-averaged flow over a period determined by a variable relaxation time on the scale of hundreds to tens of thousands of years. Measured InSAR velocities, however, are much shorter-term surface values and are subject to possible variability due, for example, to changes in basal lubrication. Differences between U _b and a measured velocity can be due, therefore, to either a short-term change in the flow behaviour or a mass imbalance due to a long-term change in the flow and/or accumulation rates. In Antarctica, no observations of significant short-term (i.e. diurnal to seasonal) variability have been reported. It is unlikely, therefore, that any discrepancy between InSAR and balance velocities will be due to short-term variability in motion.

The ability of balance velocities to reproduce both the pattern and magnitude of flow in Antarctica and Greenland has been demonstrated in two previous studies (Reference Bamber, Hardy and JoughinBamber and others, 2000a, Reference Bamber, Vaughan and Joughinb). For the Antarctic study, the comparison between U _b and the InSAR-derived velocities (U _insar) highlighted areas that showed both good agreement, and significant disagreement due to a known change in flow that took place about 130 years ago. This study demonstrated that (i) good agreement could be obtained in areas of presumed steady flow, and (ii) significant (in terms of the rms errors) differences between U _b and U _insar can be used to infer a relatively recent change in flow regime that is reflected in the surface velocity field but not in the topography, which controls the pattern of flow estimated for U _b.

InSAR velocities (U _insar)

ERS SAR data collected in late 1995 to early 1996 during the ERS tandem mission were combined to produce a velocity map of the lower reaches of PIG (Reference RignotRignot, 1998b) and TWG (Reference RignotRignot, 2001). Ascending and descending tracks of these data were combined to resolve the vector flow of the glaciers. The process involves correcting the ERS interferograms for the glacier topography, which was estimated by combining differential interferometry data with a DEM of the ice sheet (Reference Bamber and BindschadlerBamber and Bindschadler, 1997).

Velocity mapping presents two major difficulties in this region: (1) phase unwrapping and (2) velocity reference. Phase unwrapping (the process of counting phase cycles from a reference zero phase) is difficult because of the rapid rates of surface deformation of the glaciers. It was impossible to unwrap the 1992 6 day interferograms on PIG or TWG. Reference points of known velocity are required to provide an absolute reference for the InSAR measurements and to ensure that there are no residual gradients in the velocity map. In the case of PIG, the area of fast flow is relatively well confined, surrounded by numerous volcanoes and mountain peaks emerging from the surrounding ice (Fig. 1). This makes zero-referencing easy, provided that a large enough area is processed along each ERS track, i.e. several 100 km long frames. In the case of TWG, control was provided by a mountain range west of the glacier, and an isolated ice dome near the coast to the east (Fig. 1).

InSAR surface velocities for TWG were initially produced using the DEM of Antarctica, and subsequently revised using an InSAR-derived DEM of the region (from the same data). We found a difference in ice velocity of 5% between the two maps because of the better control with the InSAR DEM on areas of no motion (the Antarctic DEM was relatively poor over our areas of control, as these are isolated outcrops with large, short-wavelength relief). In areas where we have control points, we estimate that the accuracy of the velocities for the two glaciers made from ERS data acquired in 1996 is typically about 10 m a⁻¹. Inland, away from the areas of control, the errors may be higher than this. In some other areas, phase unwrapping could not be performed, resulting in gaps in the derived velocity field (shaded black in Figs 4a and 5a).

Fig. 5. (a, b) As for Figure 4, but for TWG. The maximum length of the arrows is equivalent to a velocity of 2000 m a⁻¹

Errors in balance velocities

The ice-thickness coverage includes profiles running down-flow on each glacier, as well as a single profile across TWG and several across PIG (Figs 1, 4 and 5). The thickness near the grounding lines was constrained assuming hydrostatic equilibrium for the floating ice and converting surface elevation to ice thickness. Nonetheless, significant interpolation was required to produce the regularly spaced 5 km grid produced by the BEDMAP consortium, especially for TWG (Reference Lythe and VaughanLythe and others, 2001). Accumulation data are more sparse, especially for the basin containing TWG. Errors in these two parameters scale linearly with U _b. Thus, a 10% underestimate in thickness results in a 10% increase in the estimated value of U _b. It is clear that errors in Z cannot explain the velocity differences seen at the grounding line of TWG (Fig. 7c). To investigate whether this conclusion was also supported inland, U _b and U _insar were compared along part of a RES flight-line (shown in Fig. 6), where Z is well constrained. The results of this comparison are shown in Figure 8. The very large difference (of as much as a factor 3) between U _b and U _insar from 10–50 km along track clearly cannot be explained by errors in Z. Interestingly, there is good agreement between U _b and U _insar from 100–150 km along-track, in the area outside of the fast-flowing ice. This coincides with the area of green and orange between about 103° and 105° W in Figure 6. This good agreement supports the conclusion that it is not errors in the input datasets that cause the large differences between Ub and U_insar but that a fundamental change in the flow regime is probably responsible. Given the coverage of ice-thickness data, and noting that Z is also well constrained near the grounding line (where large values of ΔU are found), it is difficult to see how errors in Z can explain the large differences seen in Figure 6.

Fig. 8. A graph of balance velocity (solid line) and InSAR-derived velocity (crosses) along the flight-line shown in Figure 6 for TWG. Also plotted is the measured ice thickness (dashed line). The track is plotted south–north, and the change in direction at about 75.5° S, 105° W takes place about 110 km along-track.

It is unlikely that errors in accumulation rate are responsible for the large differences in velocity, as this is a relatively slowly varying parameter spatially, and what is required is a large (>300%) and variable error over distances of tens of kilometres. Short-wavelength errors in surface slope could produce an offset in the position of peak flow, but TWG is in an area of relatively low surface slope (Fig. 2a) and low errors (Fig. 3b) with good coverage by satellite radar altimeter data, so this is also an unlikely explanation. It is thus difficult to explain the poor agreement in both the pattern and magnitude of flow for TWG by errors in the input datasets.

The area of negative differences for TWG cannot be explained by the variable influence of longitudinal stresses (i.e. a breakdown in the assumptions about flow). The area is spatially extensive and smoothly increasing from a minimum near the grounding line. One final possibility is that the ratio of depth-averaged to surface velocity varies dramatically but this is not supported by evidence from thermomechanical models of ice-sheet flow (Reference Bamber, Hardy and JoughinBamber and others, 2000a), which predict a ratio of 0.87–1.0 for the Greenland ice sheet. The value of this ratio is generally assumed to be close to 1.0 at the grounding line, yet, as mentioned, large values of ΔU are found here.

Errors in InSAR velocities

Velocities from descending tracks of SAR data for TWG are well constrained because they include areas of no ice motion (volcanoes, coastal peaks). The ascending tracks covering TWG have a longer interferometric baseline, which is poorer for velocity mapping, and do not have such reliable velocity control east of PIG. Control was only provided over the 50 km long sector in the lower part of the track (Fig. 1). Overall, however, and based on improvements of the velocity map compared to an earlier map, we estimate that ice velocity is measured with an accuracy of 10–20 m a⁻¹ over PIG and most of TWG, as discussed earlier.

State of balance

A short-term, historical change in accumulation rate (rather than an error in it) would cause the estimated balance flux (and consequently U _b) to be different from the present-day surface velocities. As discussed above, a very substantial, complex and implausible pattern of decrease in accumulation over much of the area covered is required, however. There is no evidence of such a change in accumulation rates from either ice cores or general circulation model simulations of the climate in Antarctica over the last century. Furthermore, the recently determined imbalance in PIG has been associated with a change in flow regime and not accumulation (Reference Shepherd, Wingham, Mansley and CorrShepherd and others, 2001). We believe, therefore, that a localized change in accumulation rate cannot fully explain the large differences seen in Figures 5, 6 and 8.

Time-invariant flow

The absolute difference in velocities for TWG (Fig. 6) indicates that the largest differences (>1000 m a⁻¹) could be reduced if the peak flow in U _b were shifted further east (upwards in Fig. 6). The discrepancy, however, covers an area of about 5000 km², and so cannot be due to the smoothing distance used in the balance-velocity calculation or to short-wavelength errors in the ice-thickness data due, for example, to an undetected, subglacial feature. The cause of the discrepancy is “broad-scale” affecting a large part of the grounding line and inland areas where ice thickness is well constrained (Figs 7 and 8).

We believe the most plausible explanation for the differences observed for TWG is a spatially extensive alteration in the flow of the glacier and surrounding ice. We propose that the glacier may have accelerated in the last few centuries, producing a mismatch between present-day observations of velocity and ones which represent flow over a longer timescale. This inference is similar to that obtained when U _b and U _insar were compared for Ice Stream C (Reference Bamber, Vaughan and JoughinBamber and others, 2000b). In this case, however, it was already known that the ice stream had recently shut down (∼130 years BP) and what was observed was the converse of the difference seen on TWG: the surface velocities from InSAR were those associated with creep flow (on the order of 50 m a⁻¹), whereas the balance velocities produced values of ∼250 m a⁻¹ (i.e. U _b was a factor five greater than U _insar).