1. Introduction
Despite several decades of intensive field investigations, we still do not know whether the Greenland and Antarctic ice sheets are increasing or decreasing in volume. The answer to this question has clear relevance to our understanding of, and ability to predict, sea-level change. The goal of NASA's Program for Arctic Regional Climate Assessment (PARCA) is to measure and understand the mass balance of the Greenland ice sheet (Abdalati and others, in press), with direct measurements of changes, over time, in ice-surface elevation as the primary objective. The program also includes “traditional” mass-balance investigations by comparing ice discharge with total snow accumulation. This is the work that we report here.
Figure 1 shows the locations of stations where ice motion is inferred from repeated global positioning system (GPS) measurements. Those along most of the western part of the ice sheet were established in 1993 and 1994 by snowmobile, and remeasured in 1995 and 1996 using a Twin Otter airplane for transportation. This part of the traverse lies within the percolation zone. We use ice velocities derived from these stations, and measurements of ice thickness along the-same route, to calculate total ice discharge for comparison with total inland snow accumulation in order to infer the rate of thickening or thinning of this part of the ice sheet. Although this is a volume-balance calculation, we use the term “mass balance” elsewhere in this paper because this is commonly encountered in glaciological literature, and volume balance can readily be converted to mass balance with knowledge of ice density.
Each traverse station comprises a jointed aluminum pole, which is the reference marker for GPS measurements, and a closely adjacent flagged bamboo pole to aid relocation of the station. Freshly planted, both poles extend about 2.5 m above the surface, and the exposed length of each pole was measured to provide an estimate of local snow accumulation. For the snowmobile traverse, one or two stations were established each work day, whereas with the Twin Otter as many as ten stations were established or remeasured in a day, with brief stops to establish a GPS receiver at the site, and a second visit later in the day lo retrieve the receiver. The GPS receivers acquired several hours of data at most stations, with a minimum occupation of about 40 min. Comparison of the repeated GPS locations, separated by a minimum of 1 year, yields estimates of surface ice velocity with magnitudes, generally, of a few tens of ma−1. Most of the GPS solutions are accurate to better than 20 cm in all coordinates, but additional errors are introduced by tilting of the marker poles and inadvertent mis-location of the GPS antenna. Here, we assume the effect of these various errors to result in a velocity error of <0.5ma-1.
Fig. 1. Greenland map showing locations of stations where ice motion is inferred from repeated GPS measurements. Contours of snow snow-accumulation rates from Reference Ohmura and ReehOhmura and Reeh (1991) in cm a water equivalent are also included. Results presented here are for the boxed area.
Airborne ice-thickness measurements were made with a coherent radar depth sounder operating at a center frequency of 150 MHz (Reference Allen, Wohletz and ChuahAllen and others, 1997; Reference Gogineni, Jezek and MooreGogineni and Others, 1998). The system transmits a chirped pulse, generating a surface acoustic wave (SAW) expander of 1.6μs dura-lion and with a peak power of 200 W, The receiver compresses the received signal in a weighted SAW processor to obtain a compressed pulse width of 60 ns. The received signal is downeconveted into baseband using a coherent detector for generating in-phase (I) and quadrature (Q) analog signals which are low-pass filtered and digitized in pairs with two eight-bit A/D converters at a rate of 18.75 MHz. Digitized signals are coherently integrated by summing 256 data samples. Integrated I and Q, signals are squared separately and summed together to determine the received power in each range gate and integrated further (incoherent integration) by summing a minimum of four consecutive samples. Data are displayed in real time to the operator and electronically recorded.
Two antenna arrays consisting of four half-wave dipoles, mounted under the left and right wings, arc used for transmission and reception. Two-way half-power beamwidth of the antenna is about 18° across the flight path and 66° along the flight path. The coherent integration of 256 pulses at a pulse-repetition frequency of 9.2 KHz reduces the along-traek antenna bandwidth lo about 17°. A sampling rate of 18.75 MHz results in the dimension of each cell being 8 m in free space or, with a refractive index of 1.78, a range-cell dimension in ice of 4.494 m. The radar thicknesses used here are estimated by computing the number of range cells between the peak signals representing the surface and the bedrock radar reflections, and multiplying ibis number by 4.494 m. During 1995 and 1996, several flight paths passed within a few hundred meters of the GRIP and GISP core sites. The radar-derived ice thicknesses agreed with the core thicknesses to within 10 m.
Estimates of snow-accumulation rate, taken from Reference Ohmura and ReehOhmura and Reeh (1991), were summed over the region corresponding to the catchment area for each of the exit gates lying between two adjacent velocity stations. These catchment areas were estimated by reconstructing the flowlines passing through all velocity stations, assuming the ice t? move in the direction of maximum regional surface slope. A surface-topography dataset compiled by Reference EkholmEkholm (1996) was used, and in what follows we shall refer to this dataset as the KMS DEM. From the KMS DEM a 1 km grid was compiled in the Universal Transverse Mercator projection system. This grid includes much of the detail of surface undulations that exist over parts of the ice sheet. Consequently, before inferring flowlines, we smoothed the data over sufficient distance to yield surface slopes likely to be those determining ice-flow direction. To accomplish this the KMS DEM was smoothed by convolution filtering using standard Gaussian filters of different sizes. The direction of maximum surface slope was calculated at each gridpoint for the original and the smoothed topographies, using data from 3 by 3 neighborhoods (3 by 3 km). For each of these cases, the direction of maximum slope was compared with the measured ice-flow direction at each of our GPS velocity stations to determine which of the smoothed topographies yielded surface slopes that best approximate the ice-flow directions.
This work resulted in selection of a 21 by 21 km Gaussian convolution filter, which removes undulations with wavelength less than about 16 km, and smooths over a distance equivalent to 10-15 times the ice thickness, which is similar to the smoothing factor used by Whillans and Reference Whillans and van der VeenVan der Veen (1993) on a West Antarctic ice stream. Flowlines were computed inland from each GPS velocity station using the smoothed topography to provide the direction of maximum slope at 5 km intervals along the flowline. Comparison between the direction of measured velocity and that of the maximum smoothed slope at the GPS stations gave a standard deviat ion of < 5°. The resulting flowlines are shown in Figure 2.
2. Analysis and Results
Within the catchment area (S), corresponding to an exit gale bounded by two of the traverse stations,
>where V is the volume of ice added each year by snow accumulation, and A is the annual accumulation rate averaged over the area S and expressed as ice thickness per year. The ice-discharge flux (F) through the gate is
where H is ice thickness, U is surface-ice velocity, and
and the integral is taken across the gate in a direction normal to ice flow.
Ice velocities were interpolated between adjacent stations, assuming linear change in speed and direction between the two measured values. The product of these interpolated values with the corresponding values of ice thickness was then integrated across the gate, with a correction factor applied to convert the gate width to its equivalent value (D) normal to the flow direction. Finally, to obtain an estimate of the ice flux through the gate, the resulting value was multiplied by a value of R derived from a model simulation of the ice sheet (Reference HuybrechtsHuybrechts, 1996) dial takes account of basal sliding and a variable temperature with depth (personal communication from P. Huybrechts, 1995). Except near Camp Century, immediately north of Rink Isbrae where the ice is moving very slowly, Huybrechts’ model shows the basal temperature to be at the melting point, and for most gates, the value of R is larger than 0.95.
Of major concern in a mass-balance calculation of this type are the errors involved in calculating the two large quantities total snow accumulation and total ice discharge— which will be compared to yield an estimate of imbalance. Consequently, we completed an error analysis which addressed these issues.
Fig. 2. Western part of the Greenland ice sheet (boxed in Fig. 1) with elevation contours in meters, showing derived flow-lines passing through velocity stations. Measured ice velocities are indicated by bold lines
2.1. Errors in measurements of ice velocity and ice thickness
Ice velocities (U) are typically 80 ma-1, and arc probably accurate to better than △U = 0.5 ma-1 in magnitude and △U/U rad in direction. Errors should be independent at each station. The errors in velocity direction introduce errors (△D) into the calculated width, normal to the ice-flow direction, of the exit gate. For most gates, △D/D is < 0.002. Measured ice thickness is typically about 1500 m, and errors probably have a systematic component (△H) of about 10 m and a random component of a similar magnitude; the effect of the random component can be disregarded, being small because of the large number of measurements made within a gate.
2.2. Errors in the assumed ratio (R) between surface and column-averaged ice velocities
The Reference HuybrechtsHuybrechts (1996) model suggests that much of the western traverse lies over ice that is not frozen to the bed, and it would be difficult to explain the observed ice velocities if this were not the case. The effect of warm basal ice and basal sliding is to increase the column-averaged velocities almost to the surface values. The model run from
which Huybrechts' R values are derived is consistent with both the present shape of the ice sheet and observed values of snow-accumulation rate, and consequently with steady-state ice velocities. Nevertheless, it is difficult to assess what error should be assigned to R; here we assumed an uncertainty of 5%.
2.3. Errors in assumed snow-accumulation rates
These are responsible for the largest error in our results, but they are very difficult to quantify. The Reference Ohmura and ReehOhmura and Reeh (1991) estimates are based on most of the data available at the time of their study, but these data are of mixed quality, were obtained over a long period and are far from uniformly distributed over the ice sheet. As new data become available, it is clear that improvements can be made. However, many of these data have not yet been published, and we have used the Ohmura and Reeh dataset because it is readily accessible to other investigators, and can be rapidly updated as new information becomes available. In the region under study, accumulation rates (A) are typically 0.4 m ice a-1, and we assume that they have errors (△A) of ± 15% that are systematic locally (10 000 km 2) but random over very large distances. Thus, averaged over very large regions (100000 km2), accumulation-rate errors probably decrease to < 10%.
2.4. Surface-topography errors
Errors in the relative elevation of neighboring gridpoints in the surface-topography dataset cause an error in the derived flowlines which affects the accuracy of our estimates of the catchment area S. For the smoothed data used here, the standard deviation between derived flow directions and those observed at the GPS velocity markers is approximately 5?????????, and we assume that this is typical along an entire flowline. Flowline directions are calculated at approximately 5 km intervals using local data from the smoothed topography dataset, and we assume that each of these calculations has a standard deviation of 5?????????. For a How band corresponding to an exit gate, this uncertainty introduces an error into the calculated width of the flow band that slowly increases with increasing distance from the exit gate. For a 350 km long flowline (typical for this study), the resulting error < △S) in catchment area is approximately 1400 km2, and this is independent of the width of the gate. For a gate between two adjacent velocity markers, S is typically 8000 km 2.
2.5. Errors in calculated thickening rates
The error (△V) in volume (V) of ice added each year by snow accumulation within a catchment area (S) is approximately:
The fractional error in V is then:
Fora gate between adjacent velocity markers △S/ S ~ 0.18, △A/A ≈ 0.15, and △V/V ≈ 0.23. The fractional error in calculated ice-discharge flux is approximately:
Using the values listed above, △H/H ≈ 0.007, △U/U^\ ≈0.006, △R/R ≈ 0.05. △D/D ≈ 0.002. and △F/F ≈ 0.05. The average thickening rate within an area S is:
with an error:
Assuming that the ice sheet is in balance, and using the values given above for an exit gate between adjacent traverse stations, S ≈ 8000km2; A ≈ 0.4ma-1; △V/V « 0.23; △F/F ≈0.05; and V = SA ≈ 3.2 x 109 m 3 ice a -1 ≈F, then △V ≈ 7.4 x 108m3 ice a-1, △F≈ 1.6 x 108 m3 ice a -1. so that:
2.6. Results
Table 1 gives values of S,A, V, D, average velocity (U) and ice thickness {Ĥ). R, F and T for each of the gates shown in Figure 2. In most cases the value of F is not identical to the product [\DUHR] because F is calculated by integrating the product [UHR] across the gate in a direction normal in ice motion. Figure 3 shows our calculated thickening rates (△T) for catchment areas corresponding to individual exit gates, with △T plotted against distance along the 2000m traverse. Actual errors are typically ±10 cm a -1, but they range from ± 6 to ±20 cm a -1,, with the largest errors for the smallest catchment areas or where is a large velocity gradient across the gate. These errors arc large, primarily because of errors in the estimated catchment area and in the accumulation rates. For larger gates, with catchment areas of 30 000 km 2, △S/S < 0.05, and part of the accumulation error becomes random, with the result that △V/V decreases to approximately 0.15, and △T ≈ ± 6 cm a-1. Consequently, we calculated thickening rates and errors for a series of gates comprising sufficient adjacent traverse stations (between 4 and 6) to represent a catchment area of about 30 000 km2. The north side of the gate was shifted, one traverse station at a time, to give the values depicted by the grey band in Figure 3.
Table. 1. Values far gates shown in Figure 2
Fig. 3. Rates of ice -thickness change (T) along part of the west side of the Greenland ice sheet. The distance, D, is measured perpendicular to ice flow from the northern end of the ice-motion traverse. Values of T for individual exit gates are shown by the vertical bars, and associated errors are between ±6 and ±20 cm a-1. The grey band depicts values of T calculated for larger gates, comprising several traverse stations, such that the catchment area is approximately· 30000km2. Here, the uncertainly in T is indicated by the thickness of the grey band, which is a smoothed version of the bar plot that is weighted according to catchment area, and consequently peaks and troughs in the two plots do not necessarily match. U and J indicate sites discussed in the text.
To reduce errors further, we calculated the mass balance for the four large regions shown in Figure 4, with areas of 80 000-100 000 km2, and with distinctive patterns of thickness change. Here, all errors are insignificant compared to those associated with uncertainty in the average snow-accumulation rates within a region, which we assume to be ±10%. This is equivalent to an error in the ice-thickness change of between ±4 and ±5 cm a -1 which represents about the highest accuracy to which we can estimate mass balance by comparing snowfall with ice discharge until we hâve better estimates of snow-accumulation rates.
Fig. 4. Four regions, A-D, with distinctive patterns of ice-thickness change shown, in cm a-1, by the circled numbers. U and J indicate sites discussed in the text, and hatching indicates the small area that appears to be rapidly thinning.
Taken together, all four of the regions shown in Figure 4 are in balance, with an average thickening rate of only 1.5 cm a -1. Fortuitously, this value is almost identical to that estimated for the entire ice sheet from a modeling study of its long-term behaviour (Huybrechts, 1994), There are two regions of thickening, B and D, separated by a region of possible thinning immediately north of Jakobshavn Isbrae.
The most stiking aspect of Figure 3 is the high variability in thickening rates for catchment areas corresponding to individual exit gales. The larger negative values are associated with very small catchment areas, and associated errors are large. In particular, the very large thinning rate of 80 cm a-1 is for the small triangular catchment area in the middle of the traverse (Fig. 4), where there is a large velocity gradient across the gate. This gate straddles the ridge separating the Upernavik drainage basin from that of the Rink Isbnr. Consequently, we calculated ice discharge through the gate by assuming zero velocity at the gate center, approximately at the crest of the ridge. With this assumption. the total volume of ice discharged was more than three times larger than the volume accumulated, but this could result from errors in our estimation of the catchment area. The neighboring gate to the south has a very large catchment area, and has a calculated thickening rate of 17 cm a-1. Taken together, the catchment areas for the two gates have a calculated thickening rate of only 5 cm a-1.
Further north, within the catchment area labelled U in Figure 3, G. Hamilton and I. M. Whillans (personal communication, 1998) measured a local thinning of more than 50 cm a-1. We also infer thinning (of 20 ± 10 cm a -1) for the single-gate catchment area containing this site (Fig. 3). Our results show both of the neighboring single-gate catchment areas to be thickening, and we originally assumed this large variability over small distances to be caused by errors in our reconstruction of the flowlines. However, the reasonable agreement with the single-point estimate of Hamilton and Whillans suggests that there may be a high degree of spatial variability in ice-thickening rates in this region. If there is, then highly localized measurements of ice-thickness change should be interpreted with great care.
Errors are significantly less for results from the approximately 30000 km2 catchment areas, depicted by the grey band in Figure 3. Here, the dominant features are the high rates of thickening close to Upernavik Isstrom (14 ± 6 cm a -1) and Jakobshavn Isbrae (10 ± 6 cm a -1). These thickening rates are equivalent to about 35% (Upernavik) and 20% (Jakobshavn} of the snow-accumulation rates averaged over the catchment areas. Evidence from recent cores (personal communication from R. C. Bales and others, 1998) indicates that Ohmura and Reeh's accumulation rates are approximately 25% too large at a site upstream of Upernavik Is-strom (U in Fig. 4) but very close to observed values at a site slightly north of the Jakobshavn drainage basin (J in Fig. 4). Thus, although the single-gate thickening rate at Upernavik in Figure 3 is probably too high, we believe our error assumptions to be reasonable, and that actual thickening rates, averaged over 30 000 km2 areas, lie within the shaded band shown in Figure 3.
Our results are broadly consistent with those from other studies. Reference WeidickWeidick (1991) presented evidence for retreat of many of the outlet glaciers in southern Greenland, from the mid-19th century until the 1950s, since when there has been a readvance. Of particular relevance to our study is the general readvance of the ice sheet south of Jakobshan Isbrae. (our region D), with continuous retreat immediately to the north (our region C). Reference ZwallyZwally (1989) and Reference Davis and HainesDavis and others (1998) interpreted satellite radar altimetry data to conclude that the ice sheet, south of 72°N, thickened between 1978 and 1986-88. Zwally’s results show significant thickening for most of this area, but recent analysis taking account of improved satellite orbits, etc., indicates average thickening rates of about 6 cm a-1, or about 25% of the earlier estimates (personal communication from H. J. Zwally, I997). Davis and others' results indicate an average thickening rate of <2 cm a-1. Their analysis of the spatial distribution of thickness change shows remarkably close agreement with our estimated thinning and thickening rate, respectively, in regions C and D (Fig. 4). Reference Krabill, Thomas, Kuiivinen and ManizadeKrabill and others (1995) compared aircraft laser data to earlier satellite Doppler position fixes on the ice sheet to infer thickening of > lm between 1980 and 1993 at locations approximately 300 km south of the region (D in Fig. 4) where we infer thickening of 6 cm a-1. There are mixed results from the so-called EGIG line, approximately along the margin between our regions C and D, where repeated leveling indicated thickening of about lm between 1959 and 1968 (Reference SeckelSeckel, 1977, but thinning by about 3.5 m between 1968 and 1992 (Reference KockKock, 1993) However, it is difficult to assess the reliability of the EGIG results. It is important to note that these various mass-balance estimates refer to different time periods, with our values indicative of trends over the past few decades.
Data collected during 1997 will permit us to complete our mass-balance study around the entire traverse shown in Figure 1. Then, during 1998 and 1999, repeat airborne surveys of all the 1993-94 laser-altimeter survey lines will reveal all areas of significant thickening or thinning over the ice sheet: and other PARCA investigations should help identify the causes of these changes. Changes detected by the airborne laser surveys will refer to the 5 year period 1993-94 to 1998-99. Longer-term trends will be measured by NASA's Geoscicnce Laser Altimeter System (GLAS), to be launched in 2001. This will accurately measure ice-sheet elevations at latitudes up to 86°, covering the entire Greenland ice sheet and most of Antarctica. GLAS data will extend the time series of Greenland aircraft measurements, and will provide a first opportunity to measure ice-thickness changes over large areas in Antarctica.
Acknowledgements
We thank M. Fahnstock for his help in remneasuring some of the traverse stations in 1995, W. Krabill and J. Sonntag for processing the GPS data, Hong Xu for assistance with the mass-balance computation, and P. Huybrerhts for providing results of his model simulations of the ice sheet. We are particularly indebted to the pilot, J, Finnbogason, and crew of the Greenlandair Twin Otter used in 1995 and 1996 to revisit the traverse stations. Their skill, patience and hard work made this project feasible. We thank R. Bindschadler and D. MacAyeal for helpful comments that significantly improved the text. This project is supported by NASA's Polar Research Program. This is Byrd Polar Research Center contribution No. 1116.
