1. Introduction
Ice sheets play a major role in studies of global climate change and its impact on mean sea level. Antarctica is the largest
reservoir of fresh water in the world, and any changes in its volume will therefore have a direct impact on the global mean sea level. To predict future changes of the ice-sheet volume, both the changes in atmospheric conditions and the dynamic response of the ice sheets require further investigation. Although our knowledge of ice physics, glacier dynamics and techniques for monitoring and modelling ice sheets has considerably improved during the last 50 years, many physical processes affecting an ice sheet still remain poorly understood. This lack of understanding is partly due to the intrinsic complexity of the behaviour of the ice but is also caused by the present lack of accurate data with global coverage, especially for the accumulation rates and ice thicknesses.
Over the last few decades, it has been established that the ice discharge from the Antarctic ice sheet is mainly controlled by outlet glaciers and ice streams. A broad description (Reference Drewry, Jordan and JankowskiDrewry and others, 1982) shows that 90% of the ice of the Antarctic ice sheet is discharged by fast-flow features which occupy only 13% of the coastline (Reference Morgan, Jacka, Akerman and ClarkeMorgan and others, 1982). Reference JoughinJoughin and others (1999), using synthetic aperture radar (SAR) techniques in the West Antarctic ice sheet region, show that tributaries can provide a spatially extensive transition between slow inland flow and rapid ice-stream flow. More recently, an estimation of balance velocity over the whole continent confirmed this result (Reference Bamber, Vaughan and JoughinBamber and others, 2000b). This evidence, which shows that ice sheets are controlled by a complex tributary-type flow, has changed our view of the dynamics of an ice sheet.
During the last 10 years, computation of ice-sheet balance velocity has become a useful tool in the study of ice-sheet dynamics. Two main computing methods have been developed: the flowline method, which is the historical method by which the balance velocity is computed (Reference Radok, Barry, Jenssen, Keen, Kiladis and McInnesRadok and others, 1982; Reference Remy and MinsterRemy and Minster, 1993; Reference Joughin, Fahnestock, Ekholm and KwokJoughin and others, 1997), and a finite-difference scheme developed by Reference Budd and WarnerBudd and Warner (1996) (Reference Bamber and HuybrechtsBamber and Huybrechts, 1996; Reference Fricker, Warner and AllisonFricker and others, 2000). The aim of this paper is to evaluate these two methods in an area where a series of high-precision global positioning system (GPS) measurements are available.
2. Computing Schemes
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:
where x is the curvilinear distance from the flow source (x = 0), Ub is the balance velocity, H is the local thickness, b(x) is the net accumulation rate and l0 is the initial distance between two flowlines (we used l0 = 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. Comparison of the Two Methods
Few comparisons have been made between the two computation schemes. Reference Bamber, Hardy and JoughinBamber and others (2000a) compared the two methods across the northeast Greenland ice stream and noted that the flowline method tends to over-concentrate the flow at the centre of the stream and overestimates the balance velocity. This comparison was made with a 20 km smoothed DEM of the Greenland ice sheet. We will investigate here the role of this smoothing distance in our area of interest and compare against in situ GPS values in section 4.
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.
4. Comparison with GPS Data
Between 1990 and 1995, Australian scientists carried out a series of glaciological traverses around the LGB. This traverse route, between the Mawson and Davis stations, approximately followed the 2500m surface elevation contour, and data were collected at 73 ice-movement stations spaced approximately 30 km apart. Recently, Reference Manson, Coleman, Morgan and KingManson and others (2000) have processed the GPS data collected at these 73 stations along the LGB traverse route to produce precise surface velocities. Before the GPS velocities can be compared with the balance-velocity results, the GPS surface values must be converted to depth-averaged velocities. The ratio between surface speed and depth-averaged speed depends on ice rheology and the physical conditions in the ice column and at the base ofthe ice layer. A ratio of 0.87 has been applied between depth-averaged and surface velocity values, as defined in Reference Budd and WarnerBudd and Warner (1996).
Figure 4 presents the comparison between the 73 GPS velocity measurements (black line) and the 2161 computed velocity points along the LGB (grey lines) using the 5 km DEM (Reference Remy, Shaeffer and LegresyRemy and others, 1999) with L =100 km. Figure 4 shows that both computed balance velocities are in good agreement with the GPS-derived estimates. However, a closer look at Figure 4 also reveals that significant differences still exist, especially in the high-flow glacier regions, such as Mellor and Lambert Glaciers. In order to assess the capability of both methods to reproduce the intensity along the LGB traverse route, we have defined six separate zones (see Fig. 4). In each of these zones, the mean relative difference between modelled and observed velocity was computed as a function of the smoothing length scale (see Fig. 5). We adopted the convention that all percentages shown in parentheses below correspond to the FL method.
Fig. 4. Balance velocities around the LGB traverse route. The solid black line is the GPS velocity, the dashed grey line computed velocity and the solid grey line the FL velocity. The numbers along the black line are the LGB stations, starting western side of the basin and finishing on the eastern side.
Fig. 5. The relative difference (%) between computed and GPS velocities for the whole LGB traverse profile and for each of the six zones defined in Figure 4. The solid black line is the FL velocity, and the dashed black line the BW velocity. The grey lines are the corresponding relative differences when sliding is considered (ratio = 1 for Vgps > 30 m a–1).
For the whole profile, the mean computed balance velocity is underestimated by about 9% (15%) compared with the mean GPS velocity. In low-motion regions (zones 1, 4 and 5) we can see (Fig. 4 and 5) that the balance velocities are quite well modelled: the computed velocities in the western part of the LGB (zone 1) are higher by ∼10% (2– 3%) when compared with the GPS values, and in the eastern part (zones 4 and 5) the values are generally smaller by ∼10% (12%) and ∼5% (10%) respectively. These percentages are within the input data uncertainties and methodology, but a small imbalance between these two sides of the LGB could also explain these differences. The high-velocity regions (zones 2,3 and 6) are found to have higher computed balance velocities of about 25% (15%) for Lambert Glacier, near-identical velocities 3% (–3%) for Mellor Glacier and smaller velocities of about 18% (–20%) for the eastern part. A part of this difference in velocities within the higher-velocity regions could be due to ice-sheet sliding. Indeed, when sliding is assumed for velocities >30ma–1 (i.e. for these velocities, the ratio between depth-averaged and surface velocity is taken to be 1instead of 0.87), the relative difference between computed and measured velocities decreases to 10% (1%) for Lambert Glacier (but increases for zones 2 and 6). The spatial sampling of the LGB stations (i.e. distances of about 30 km between stations) could also miss some high-velocity, smaller-scale streams, which would lead to an under-estimation of the mean GPS velocity.
We can see from these comparisons that the methodology adopted for the balance-velocity computation accounts for 5– 10% of the error in the final velocity. We review other possible sources of uncertainty in the next section.
5. Other Sources of Uncertainty
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.
Conclusion
Comparison between computed balance velocities and GPS measurements, observed along the Australian LGB traverse route, shows that actual balance-velocity computing schemes are able to recover the spatial pattern of the velocities and to model their average intensity to within 5– 25%. The major sources of error in this region, where ice thickness values have been collected and hence are assumed well known, are found to be the uncertainties present in the accumulation rate, the adopted computational scheme (either the FL or BW method) and the smoothing length scale (L) applied to the topography. With improvement in the resolution of Antarctic DEMs likely from satellite missions such as IceSat, and an increasing refinement in the accumulation-rate and ice-thickness data, the dominant sources of these errors can be removed. Despite improvement in input data sources, the role of the filter used in the balance-velocity calculations still needs to be further investigated, especially aspects such as the type of filter used, the optimum smoothing distance, etc.
In this study, the BW method was found to be very robust, and with high-resolution DEMs (1km grid size) the method was able to capture the high-flow regions of the glacier streams. The FL method was very sensitive to the adopted smoothing-length scale, but when the topography was sufficiently smoothed (at longer spatial scales), the resulting balance velocities were equivalent to the BW values.
Acknowledgements
The work by R. Hurd and R. Coleman was supported by grants from the Australian Research Council and the Australian Antarctic Science Grants Scheme. R. Hurd was supported by an Australian Postgraduate Award held at the University of Tasmania. The authors are grateful to the scientific editor and the two anonymous referees for their valuable comments which led to great improvements in this paper.
