Assumptions

A classic way of estimating the drainage pattern in Antarctica is to compute the so-called “balance velocity”. The balance velocity corresponds to the depth-averaged velocity needed at any point to maintain the ice sheet in a state of balance. Two assumptions are required to make the computation.The first one presumes that the ice sheet is in a state of balance, such that the outward flow distribution matches the net accumulation without considering any time variation in the thickness. The second assumption considers that, in a first approximation, the surface slope corresponds to the dominant direction of the ice flow. Under these assumptions, the mass-balance equation of an ice sheet can be written:

(1)

(2)

where x is the curvilinear distance from the flow source (x = 0), U_b is the balance velocity, H is the local thickness, b(x) is the net accumulation rate and l₀ is the initial distance between two flowlines (we used l₀ = 5 km). Hence, with a given distribution over the ice sheet of net accumulation rate and ice thickness, the balance velocity can be estimated.

Flowline and finite-difference schemes

The flowline scheme (now called FL) is based on the upstream calculation of the surface between two distinct flowlines, originally lying on a contour on both sides of a point (Reference Remy and MinsterRemy and Minster, 1993). The upstream paths of the two flowlines are then computed at 5 km intervals along the flowlines, with a local slope calculated over a 10 km length scale. In this way, each point along the flowline is computed separately from the others and is free to evolve through the specified grid resolution.

The finite-difference scheme is based on Reference Fricker, Warner and AllisonFricker and others (2000), which is an improvement of the Reference Budd and WarnerBudd and Warner (1996) method (now called BW). This method is based on a two-dimensional, finite-difference scheme, which distributes the accumulation to the eight neighbouring cells according to the slope direction. The computation is carried out for each cell, previously sorted in decreasing order, so that each cell of the grid receives the flux of all of the previous cells.

Input dataset

Computation of the balance velocities was performed in the Lambert Glacier basin (LGB) region, using the two methodologies discussed above. Two different 5 km digital elevation models (DEMs) of the Antarctic ice sheet, derived from the altimetric data of the European Remote-sensing Satellite (ERS-1) (Reference Liu, Jezek and LiLiu and others, 1999; Reference Remy and MinsterRemy and others, 1999), were chosen to assess the variability due to topography. The BEDMAP thickness compilation (Reference Lythe and VaughanandLythe and others, 2001) and a recent compilation of accumulation rates (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999) were used as the other data-input sources.

3.1. The role of the smoothing scale

The question of what spatial smoothing should be applied to various data sources when we undertake any flow computations is problematic. From a methodological point of view, under certain smoothing length scales, neither method isable to produce a coherent picture of the flow because each is sensitive to small-scale undulations in the topography, i.e. the presence of numerous isolated depressions. From a glaciological point of view, as the effect of longitudinal stresses is not included under the assumption of downslope flow, a classical smoothing distance of 10–20 times the ice thickness is applied to the topography in order to remove the effect of longitudinal stresses (Reference PatersonPaterson, 1994). However, this length scale is highly dependent on the ice rheology and flow conditions.

The approach adopted here is to study the sensitivity of both schemes to a chosen smoothing length scale applied to the topography. Heights in the Lambert Glacier drainage basin region were extracted from the available DEMs and smoothed with a Gaussian moving-average box of varying size L. The DEM was smoothed from 5 to 25 times the grid-cell size (i.e. from 25 to 125 km), leading to 21 different versions of the DEM. We applied the BW method on each of these versions to produce 21 different velocity maps. The next step was to digitize coordinates (2161 values), based on the LGB traverse (see section 4), to produce a 2000 km profile around the LGB with approximately a point every 1km. At each of these points and for each smoothed version of the DEM, the balance velocities were computed using the FL method, and the BW velocities were interpolated at these points from the corresponding velocity map.

Along this profile we have defined three different zones: zone 1 corresponds to the first 925 km of the profile (the western part of the LGB); zone 2 corresponds to the glacier inflow region, including Mellor and Lambert Glaciers; and zone 3 corresponds to the other end of the profile (the eastern part of the LGB; see Fig. 2a). For each smoothing distance and for each zone, the mean differences and the rms between balance velocities obtained with the two methods (UbFL – UbBW) were computed. Moreover, in order to see the convergence of each method, the mean difference and rms were also computed between the velocity profile obtained with a smoothed version of the DEM and the velocity profile obtained from the previous smoothed version (i.e. DEM(i)–DEM(i–1), where i = 1, 21). One can see in Figure 1 that, in terms of the mean difference and rms, the FL method (greyline) and the BW method (dashed grey line) converge towards a stable solution as the smoothing distance increases. The BW scheme is, however, found to be much less sensitive to the smoothing distance than is the FL scheme. Indeed, in low-velocity regions (i.e. zones 1and 3), until a smoothing distance ofat least 50 km is reached, the

Fig. 1.

Mean velocity difference and its rms variation. The solid black line is the difference between the FL balance velocity BW balance velocity. The grey lines (full = FL, dashed = BW) are the differences between computed velocities associated the smoothing distance L and velocities associated with the smoothing distance (L – 5 km), and show the convergence methods. The three upper plots represent the rms of these differences computed for the three zones, and the three lower mean difference (in m a^–1). The star in each box shows the corresponding quantity computed using a 1km DEM.

Fig. 2.

The three curves show the computed balance velocity for the digitized profile along the LGB traverse route for smoothing or grid sizes. Dashed grey lines represent the BW computed velocity, and solid black lines are used for the (a) and (b) are computed using a DEM grid size of 5 km, and (c) uses a DEM grid size of 1km.

FL method has much higher variability than the BW method. In high-velocity regions (i.e. zone 2), if one uses the FL method, the DEM needs to be smoothed over distances of 100 km or more to reproduce a variability comparable to the BW method. In Figure 2 we have represented the computed balance velocities along the whole profile for two selected smoothing length scales (L = 50 and L = 100 km). We can see in this figure that the FL scheme produces, at short smoothing length scales and in particular glacier regions, narrow and high-intensity peaks. These peaks are not computing artefacts, they denote the tendency of the FL scheme to over-concentrate the flow at the centre of topographic depressions. We defined the quantity

where x is the curvilinear distance from the flow source (x = 0) and l is the distance between two adjacent flowlines. This quantity is the ratio between the drainage surface of the ice sheet and the surface drained by the equivalent parallel flow and can be seen as the relative state ofconvergence or divergence of the topography. If we map this quantity in the LGB region, using the 1km version of the Reference Liu, Jezek and LiLiu andothers (1999) DEM (filtered at a 50 km smoothing scale) and using the FL scheme, we can observe that this scheme produces a very coherent pattern of highly convergent, narrow channel flow that obviously represents a signature of the topography (see Fig. 3). This signature cannot be directly interpreted as

Fig. 3.

The map represents the F index computed in the region using the Reference Liu, Jezek and LiLiu and others (1999) 1km DEM with L 50 km.

LGB a velocity map but it is most probably linked to actual or = past flow conditions.

3.2. The role of the grid size

As the BW scheme is dependent on the chosen grid size, the two methods are not exactly comparable. Indeed the grid size in the BW method plays the role of an additional smoothing parameter compared with the FL method, which is less sensitive to it. We have computed the two balance velocities along the LGB traverse profile using the 1km DEM. Figure 2c shows that, with this grid size, the BW velocity contains more high-frequency information and is now much closer to the FL velocity. This is confirmed by the mean difference and rms value (see star in Fig. 1). Thus, comparison between the two methods shows that the FL scheme is very sensitive to the smoothing filter applied to the DEM, and indeed needs a large amount of smoothing to be comparable to the BW scheme. On the other hand, the BW scheme is sensitive to the size of the gridcell.

Ice thickness

In the mass-balance equation (Equation (1)), the relative error in balance velocity due to an error in the local ice-thickness value is simply equal to the relative error in the ice thickness. This error in relative ice thickness can reach a few tens of per cent in regions where data are sparse. In the LGB region, many ice-thickness measurements have been collected, and these data are included in the BEDMAP compilation. Therefore along the LGB traverse, the relative ice-thickness error is estimated to be only a few per cent, leading to an equivalent few per cent error in the computed balance velocity.

Accumulation

The accumulation is an integrated quantity, and hence a random fluctuation in the accumulation rate will not have a great impact on the balance velocity. A fluctuation of about 10% in the accumulation rate randomly distributed on the whole ice sheet and decorrelated at the scale of the grid size (5 km) leads to relative variation in the computed velocity of only a few per cent. On the other hand, a bias in the accumulation at a regional scale will have a strong effect in terms of the computed velocity magnitude. The use of two different accumulation datasets (e.g. Reference Radok, Brown, Jenssen, Smith and BuddRadok and others, 1986; Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999), which show strong regional-scale differences, can lead to differences of >50% in the computed magnitude of the velocity.

At the LGB stations, we have performed a series of balance-velocity computations using the four “best”-performing accumulation datasets identified by Reference Fricker, Warner and AllisonFricker and others (2000) in this region: (i) the Australian Antarctic CRC dataset(Reference Budd, Reid and MintyBudd and Smith, 1982); (ii) the Global Atmospheric Assimilation and Prediction Scheme (GASP) (Reference Budd and SmithBudd and others, 1995); (iii) the Commonwealth Scientific and Industrial Research Organisation (CSIRO) dataset (Reference Smith, Budd and ReidSmith and others, 1998); and (iv) the Vaughan dataset (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999). The CRC and Vaughan datasets are based on field observations, whilst the GASP and CSIRO datasets were determined from atmospheric models. The results of these comparisons show that using accumulation values from an atmospheric model, as opposed to an observational dataset, produces velocity differences of around 30%, with significantly higher percentage differences at a number of sites (see Reference HurdHurd, 2002). In general, the choice of accumulation input datasets caused the computed velocities to vary from the GPS-derived values by up to 55%.

Surface topography

Two different altimeter-derived DEMs (Reference Liu, Jezek and LiLiu and others, 1999; Reference Remy, Shaeffer and LegresyRemy and others, 1999) have been used to test balance-velocity results in the LGB region.The computation of the 1km velocities along the LGB traverse using these two different DEMs shows, at equal smoothing, near-identical results in terms of velocity magnitude but exhibits small distortions in spatial geometry. Cross-correlation of the computed velocities shows generally very high correlation (40.92), but also indicates some lags of a few km, mostly in regions of high convergence–divergence. These lags are due to geometrical differences in the surface topographies of the two DEMs. These distortions are of the same magnitude as those that can be seen in Figure 4 between the GPS and the computed balance velocities (see, e.g., LGB stations 6 and 50). Indeed, one can see that the maximum (minimum) computed velocity occurs a few km before (after) the maximum (minimum) GPS values. Hence these lags between observed and modelled velocity can probably be explained by the fact that the DEMs used here are not exactly identical as they do not match the true surface of the ice sheet. Additionally, the smoothing applied to the topography causes modification of the shape of the surface contours.