Introduction

The large ablation zone of West Greenland (Fig. 1) is a potentially significant contributor to future sea-level rise because CO₂-induced climatic warming may increase both the rate and area of ablation. In this paper we investigate the extent to which satellite radar-altimeter data can be used for reliable measurement of average changes in the average height of the ice-sheet surface within the ablation zone. The data are from the Geosat Geodetic Mission (GM), which extended from spring 1985 through summer 1986. Altimetry offers the advantage of comprehensive coverage from the southern tip of Greenland to 72°N, which is the maximum northern extent of the orbits of altimeter-carrying satellites at present. The measurement noise levels are high and spatially variable within the ablation zone because of relatively steep surface slopes and rough topography.

Semivariogram methods are used to estimate noise in the altimeter measurements as a function of position. The semivariogram, which is a basic component of the theory of regionalized variables (Matheron, Reference Matheron1963), is analogous to the autocorrelation function as used in time-series analysis (e.g. Olea, 1974, Reference Olea1977; Curran, 1988). The noise estimates include measurement noise from all sources, and do not depend upon a priori knowledge of instrument calibration, radial orbit error, or errors from other sources (Hillger and Vonder Haar, Reference Hillger and Vonder1988).

Seasonal mean changes in the surface height are determined from the height differences at orbit cross-over points. A cross-over point is a location where an orbit ascending in latitude is later crossed by an orbit descending in latitude, or the reverse, so two measurements of the surface height — separated by a time interval — are obtained at approximately the same location.

The cross-over analysis is a variation of the method developed by Zwally and others (1989), but differs in that semivariograms are used to estimate the measurement noise as a function of position. This approach makes it possible to suppress the influence of cross-over differences from high-noise areas by weighting each cross-over difference in proportion to the inverse square of the noise level in the immediate vicinity of the cross-over point. The noise-induced errors are propagated through the cross-over analysis using standard methods (Reference Moffitt and BouchardMoffitt and Bouchard, 1965, p. 163–68).

Preliminary Corrections and Data Selection

Altimeter data over the Greenland ice sheet from the Geosat GM were reduced using the Navy precision orbits, and were corrected for atmospheric effects and solid Earth tides as described for Seasat altimetry by Zwally and others (Reference Zwally, Bindschadler, Brenner, Martin and Thomas1983) and for tracking errors using the method described by Martin and others (Reference Martin, Zwally, Brenner and Bindschadler1983).

Significant error is caused by slope effects (Reference Brenner, Bindschadler, Thomas and ZwallyBrenner and others, 1983). The absolute accuracy of the altimeter-derived elevations would be increased if corrections for slope-induced errors were applied, but this was not done because of the possibility that additional noise would be introduced. Here we attempt to maximize the accuracy of measured mean changes in the surface height, not the absolute accuracy of the elevations. Toward this end we assume there is constant correlation between the “altimeter surface” and the “real surface” at successive times, and the changes will be measured more accurately if the data are not corrected for slope-induced errors.

Data were selected by applying elevation criteria to the Geosat GM data from the western Greenland ice sheet. The upper limit of the ablation zone was considered to be the firn line plotted as a function of latitude by Benson (Reference Benson1962, p. 74). South of 66°N, the firn line was assumed to lie at the elevation of Benson's southernmost data point (1700 m). (The more recent data of Arnbach (Reference Ambach1963) and Thomsen (Reference Thomsen1984) suggest that Benson's equilibrium line may be somewhat high.) The lower elevation limit for acceptable data was defined, also as a function of latitude, at sufficient distance up-glacier from the margin to exclude wave forms back-scattered from coastal mountains and nunataks. Application of these criteria resulted in the data set shown in Figure 1.

Fig. 1. Map showing positions of re-tracked wave forms from the Geosat GM along orbit ground tracks throughout the West Greenland ablation zone.

Semivariograms of Altimeter-Derived Elevations

The surface of an ice sheet is spatially continuous but includes a stochastic component; it is thus a regionalized variable. The discrete altimeter-derived surface heights z _i throughout the localized neighborhood around a grid point can be written:

where m is the smoothly-varying drift throughout the neighborhood of the grid point, r _i is a residual which includes both the stochastic component of the surface height and the error of the. measurement, and ϕ, λ are north latitude and east longitude, respectively.

The drift m is defined to be the biquadratic surface which is a best fit, in an (unweighted) least-square sense, to the z _i throughout the neighborhood of the grid point. See Zwally and others (1990) for a detailed description of the method used to fit m to the z _i. Removal of the drift results in residuals r _i which are expressed, henceforth, as functions of orthogonal coordinates x,y on a tangent polar Stereographic projection:

where X _i = (x,y)_i. Details of the transformation from ϕ,λ to x,y have been given by Zwally and others (Reference Zwally, Brenner, Major, Martin and Bindschadler1990) and Parkinson and others (Reference Parkinson, Comiso, Zwally, Cavalieri, Gloersen and Campbell1987, p. 248–19).

The residuals, which are also a regionalized variable, have zero mean by definition of the least-squares fit, and zero linear trend because the linear component of the drift is included in m(x,y).

The semivariance of the elevation residuals at lag h _P is:

where h _p = ap, p = 1,2,…, and a is the distance between successive altimeter wave forms along the sub-satellite tracks (662 m). A tolerance δ = a/2 is defined such that all residual pairs separated by h _p ± δ have lag h _p. The tolerance δ is necessary because the measurements are not always separated by integer multiples of a, due to multiple non-parallel sub-satellite tracks crossing the neighborhood of a given grid point, and also due to gaps in the data caused by the roughness of the surface and local changes in surface slope.

Note that, although Equation (3) applies to residuals distributed arbitrarily in two horizontal dimensions, it is identical in form to the expression for average semivariance applicable to residuals at equally spaced intervals along a line (Olea, Reference Olea1974).

An experimental semivariogram is computed from Equation (3) for the circular neighborhood around each point of a 10 km grid throughout the West Greenland ablation zone, to 72°N. The neighborhoods have a minimum radius of 15 km. Figure 2 shows a respresentative semivariogram, computed for the neighborhood of a grid point located at 69.587°N, 49.699°W.

Two features of semivariograms, the range and the sill, are labelled in Figure 2. Discussions of these features are given by, for example, McBratney and Webster (Reference McBratney and Webster1986). The range is the finite lag within which the maximum mean-square variation is encountered. Elevation residuals separated by lags less than the range are spatially correlated. Those separated by lags greater than the range are spatially independent. The sill represents the a priori variance of the residuals.

Altimetry Noise Levels from Semivariograms

The average measurement error (noise) throughout the neighborhood of each grid point is determined by extrapolating the corresponding semivariogram to zero lag. The result is one-half the mean-squared difference between repeated altimeter measurements of elevation at the same geographic location. That is, the square root of the semivariance at zero lag is the standard error of the measurements throughout the neighborhood used to compute the semivariogram (Burrough, Reference Burrough1986, p. 155; Hillger and Vonder Haar, 1988).

The zero-lag intercept is located by fitting a polynomial to the semivariances for lags less than 4 km (see Fig. 2). A maximum lag of 4 km is used because it is found to be a conservative estimate of the range of most semivariograms. The polynomial is employed only for extrapolation to zero lag, and does not otherwise serve as a model of the experimental semivariogram. (Only a limited class of functions is acceptable for modeling experimental semivariograms; see e.g. McBratney and Webster, 1986.)

The polynomial is fitted using the weighted least-squares method employed by Hillger and Vonder Haar (1988). Each semivariance is weighted in proportion to the number of elevation pairs used in its computation (i.e. those computed with a larger number of pairs are weighted more heavily), and in inverse proportion to lag (i.e. semivariances near zero lag are also weighted more heavily).

Fig. 2. Representative semivariogram of elevation residuals from the ablation zone, showing approximate range and sill, polynomial (solid line) used to extrapolate to zero lag, and mean-square measurement noise throughout the neighborhood of the grid point at the stated latitude and longitude.

The criteria for the polynomial coefficients employed by Hillger and Vonder Haar (1988) are also used. (1) The coefficient of the linear term (lag to the first power) is required to be zero. The first derivative of the polynomial is thus zero at zero lag. (2) The coefficient of the “lag-squared” term is required to be positive. The polynomial thus increases with lag for small lags. This criteria also ensures that the minimum of the fitted polynomial is at the zero-lag intercept. (3) The zero-lag intercept must be positive, since the noise level is the square root of the intercept.

Fourth-order polynomials are used unless the above criteria cannot be met. In such cases the fit is attempted using a third-order polynomial. Second-order polynomials were found to yield poor fits and were not used. If the above criteria cannot be satisfied with a third-order polynomial, a larger neighborhood containing more elevation measurements is defined, a new semivariogram is computed (after removing the drift from the enlarged neighborhood), and the polynomial-fitting procedure is repeated. In general, it was found that semivariograms computed from larger numbers of residual pairs were more likely to be usable for determination of measurement noise by extrapolation to zero lag.

A map showing altimeter-measurement noise levels throughout the West Greenland ablation zone, to 72°N, is shown in Figure 3. The map is contoured from the 10 km grid. Each grid point value is the average noise level throughout the circular neighborhood (with a minimum radius of 15 km, centered on the grid point) used to compute the corresponding semivariogram. The noise levels represent the effects of all factors contributing to scatter in repeated measurements at the same geographic location,

Fig. 3. Altimetry noise levels throughout the West Greenland ablation zone, extrapolated from semivariograms. A semivariogram was computed from data in the neighborhood of each point of a 10 km grid throughout the region shown in Figure 1. The distance scale is stretched west to east (20 km per grid interval, versus 100 km per grid interval south to north).

including (i) RMS orbit error, (ii) systematic bias between ascending and descending orbits (discussed in the next section), and (iii) changes in the elevation of the ice surface during the Geosat GM. The semivariogram approach described above thus yields the desired result: altimeter-measurement noise levels as a function of position.

Elevation Changes from Cross-Over Differences

Changes in the height of an ice-sheet surface can be measured, in principle, by taking the difference between successive measurements at orbit cross-over points. In practice, however, some cross-over differences are unrealistically large.

This problem is ameliorated by editing (discarding) the cross-over differences that are larger than a defined maximum. Preliminary analyses with a relatively tight edit level (10 m) suggested that seasonal mean changes in the surface height in the ablation zone are likely to be on the order of 1-2 m. Superimposed on this are measurement noise levels in the 6-16 m range over large areas, as shown in Figure 3. The elevation at a cross-over point is obtained by linearly interpolating between the closest elevations on either side, so the standard error of a cross-over difference is (2)^½ e, where e is the standard error (noise level) for measurements in the vicinity of the cross-over point. If typical noise levels are about 11 or 12 m, and if seasonal changes in the surface height can be about 2 m, then this could result in apparent cross-over differences (signal plus noise) of roughly 15–19 m. 20 m was thus chosen as a reasonable cut-off. The 20 m edit criterion results in 7769 cross-overs being selected for use out of a total of 9993.

The measurement noise levels are known as a function of position (Fig. 3) and are calculated on a 10 km grid. The influence of cross-over differences from areas characterized by exceptionally high noise levels is suppressed by weighting each cross-over difference in proportion to the inverse square of the noise level at the nearest grid point.

Apparent height changes determined from cross-over differences include the effects of radial orbit error. In the case of Geosat data, systematic orbit bias is a significant component of the radial orbit error. The orbit bias is the difference between the surface height measured during an ascending orbit, and the height measured during a descending orbit at the same location at a different time, if the true surface height does not change during the time interval and if other sources of error do not exist. The radial orbit error is approximately 0.8 m. Of this, the average ascending-descending orbit bias in the Geosat data over Greenland is −0.46 ± 0.01 m (Reference ZwallyZwally and others, 1989). The ascending-descending orbit bias is eliminated by averaging the cross-over differences between seasons as described below.

Throughout an initial time interval Δt ₁ which is 90 d following the start of data from the Geosat GM, the measurements made during the ascending orbits are selected. The start of the data is day 91 of 1985, i.e. the beginning of April. This initial 90 d interval is defined as spring 1985. Throughout a second time interval Δt ₂, which is the second 90 d (summer 1985), the measurements made during the descending orbits are selected. Throughout the region the “summer descending orbits” cross the “spring ascending orbits” at points i = 1, 2 …, n. (These are assumed to be the n points having usable data, after editing.) At the ith cross-over point the change in the surface height, including measurement error, is Δh _i and the orbit bias is b _i. The weighted-average height change measured by these n cross-overs throughout the region is:

where:

and:

where e_i , the noise level in the neighborhood of the ith cross-over point, is taken from a digitized version of the map shown in Figure 3. Equation (4) is an unbiased average (the weights sum to 1).

Similarly, the height measurements made during descending orbits during the first 90 d are subtracted from those made during ascending orbits during the second 90 d, at the orbit cross-over points. There are now n′ cross-over points ([Inline 1]) which occur at different locations. The weighted-average height change measured by these n′ cross-overs throughout the region is:

where:

and:

In Equation (5) b has the same sign as in Equation (4), because measurements made during descending orbits are subtracted from those made during ascending orbits in each case. However, Δh is of opposite sign in Equations (5) and (4) because the surface height at an earlier time is subtracted from the height at a later time in Equation (5), while the reverse occurs in Equation (4).

The average orbit bias in Equation (5) is the same as the average orbit bias in Equation (4), if n and n′ are sufficiently large and if the orbit cross-overs provide representative coverage in each case. These criteria are assumed to be met. Multiplying Equation (4) by -1 and taking the average:

In Equation (6), the orbit bias has cancelled and the change in the surface height remains.

The propagation of the measurement noise, as determined from the semivariograms, through the cross-over computations given by Equations (4) through (6), is computed using standard methods (Reference Moffitt and BouchardMoffitt and Bouchard, 1965, p. 163–68). The standard errors corresponding to Equations (4) and (5), respectively, are:

and the standard error of the mean elevation change given by Equation (6) is:

The seasonal mean changes in the surface height relative to spring 1985 are obtained, finally, by referencing all changes between successive seasons back to spring 1985 and averaging again. For example, the change between spring 1985 and summer 1986 is computed as follows. (I) Spring 1985 is “crossed” with summer 1986. (2) Summer 1985 is “crossed” with summer 1986, and the result is corrected for the change between spring 1985 and summer 1985. (3) Fall 1985 is “crossed” with summer 1986, and the result is corrected for the change between spring 1985 and the mean surface of fall 1985, etc. The changes from spring 1985 and successive seasons to summer 1986 are thus with respect to a common reference (spring 1985). These changes are then averaged, with each change again weighted, in proportion to the inverse square of its standard error. The errors of the weighted averages are computed in a manner analogous to Equation (7). All of the relevant cross-over differences between successive seasons are thus included in the computation of each seasonal change, and the standard errors of the means are reduced somewhat further by weighted averaging.

Seasonal Elevation Changes in West Greenland

Figure 4 shows seasonal mean changes in the surface height of the West Greenland ablation zone (to 72°N), between spring 1985 and summer 1986, computed using the method described in the previous section. The data begin on day 91 of 1985, which is at the beginning of April. Spring 1985, as labelled in Figure 4, is the subsequent 90 d period (April, May, June). Summer is the 90 d period from the beginning of July through the end of September, fall is the 90 d period from the beginning of October through the end of December, etc.

The mean surface height decreased by 1.5 ± 0.6 m between spring and summer 1985. Relative to spring 1985, the mean surface height during fall 1985, winter 1986, and spring 1986 was higher by 1.3 ± 0.6 m, 0.5 ± 0.5 m, and 1.7 ± 0.4 m, respectively. During summer 1986, the mean surface height was lower by 0.98 ± 0.3 m than during spring 1985.

Comparison to Field Studies and Other Work

According to R.J. Braithwaite (personal communication), over most of West Greenland the 1984–85 winter was

Fig. 4. Mean seasonal changes in the surface height of the West Greenland ablation zone to 72°N, from spring 1985 through summer 1986, relative to spring 1985 (datum). The error bar for each season is the standard error of the mean.

relatively warm with little snow accumulation, and the 1985 summer had relatively heavy melt. Subsequently the 1985–86 winter had heavy snowfall. The average increase of the surface height of the ablation zone between spring 1985 and spring 1986 shown in Figure 4 may thus be a result, wholly or in part, of above-normal snowfall during winter and spring 1985–86.

A summary of mass-balance measurements from 1980-81 through 1983–84 along a line of stakes extending up the center line of Qamanârssûp sermia (a marginal glacier in south-west Greenland at 64.5°N) from the terminus to the equilibrium line is given by Braithwaite (Reference Braithwaite1986). When plotted, the data in his Table I (1986, p. 51) show a trend toward less negative mass balance at all stakes. However, whether this glacier was thinning, thickening, or in equilibrium prior to the study is apparently unknown, so the ablation zone may have either thickened or thinned at a decreasing rate during the measurement period.

According to Thomsen and others (1986), an area of the West Greenland ablation zone near the margin north-east of Jakobshavn, centered at approximately 69.43°N, 50.17°W (Fig. 1), thinned by an average of 14 m between 1959 and 1982. This corresponds to a mean thinning rate of 0.61 m year^×1. The area is about 22 km south-west-north-east, parallel to the margin, and extends about 6 km into the ablation zone perpendicular to the margin. Our methods cannot be applied to this area for direct comparison because it is too small to yield enough orbit cross-over points for measurement of a meaningful average elevation change. Thomsen and others (Reference Thomsen, Thorning and Braithwaite1986) also pointed out that the margin is advancing south of Jakobshavns Isfjord, at 68.75°N. This latter observation is more consistent with the results of this study, and suggests that the thinning trend observed north-east of Jakobshavn may not be representative of the entire West Greenland ablation zone.

Zwally and others (1989) found that the average surface elevation of the Greenland ice sheet south of 72°N increased during the 1985-86 time frame of the Geosat GM, and also during the 1978-85 Seasat-to-Geosat GM time interval. South of 65°N, the average surface elevation increased during the 1975–78 GEOS-3-to-Seasat time interval as well. Increasing surface elevations were measured within all elevation intervals, including the 700–1200 m interval, which generally coincides with the ablation zone. These findings are consistent with our measurement of an increase in the mean surface height of the West Greenland ablation zone between spring 1985 and spring 1986, and also between summer 1985 and summer 1986.

The error limits determined in this study are larger than those determined by Zwally (1989) for the 700-1200 m elevation interval. The semivariograms used to determine noise levels in this study were computed from all of the data in the neighborhood of each grid point (except that a 3ϕ edit was applied when fitting the biquadratic surfaces used to remove the local drift — i.e. elevations differing from the surface by more than three standard deviations were edited during the least-squares fitting procedure). Zwally (Reference Zwally, Brenner, Major, Bindschadler and Marsh1989) determined error limits by computing the mean and standard deviation for all cross-over differences throughout the 700-1200 m elevation interval. The cross-over differences greater than either 3ϕ or 10 m from the mean (whichever was less) were then edited, and the standard deviation of the remaining cross-over differences was taken as the standard error. That is a tighter editing procedure than employed here.

Summary and Conclusions

The results of this study show that semivariogram methods are suitable for estimation of satellite radar-altimeter noise levels as a function of position over the polar ice sheets. This approach permits computation of the propagation of the noise-induced errors associated with individual cross-over differences through the cross-over analyses used to determine time-dependent mean changes in the surface height. In addition, this approach permits minimization of the influence of cross-over differences from areas characterized by exceptionally high noise levels by weighting each cross-over difference in proportion to the inverse square of the noise level in the neighborhood of the cross-over point.

Relative to spring 1985, the mean surface height of the West Greenland ablation zone to 72°N was found to be lower by 1.5 ± 0.6 m in summer 1985, higher by 1.3 ± 0.6 m in fall 1985, higher by 0.5 ± 0.5 m in winter 1986 (i.e. starting at the beginning of January), higher by 1.7 ± 0.4 m in spring 1986, and lower by 0.98 ± 0.3 m in summer 1986 (see Fig. 4), using data from the Geosat GM.

These results suggest, in particular, that the method will be suitable for measurement of mean changes in the surface height near the margins of the polar ice sheets, where altimeter noise levels are high and variable due to steep surface slopes and rough topography, and that it will be possible to measure longer-term elevation changes that may be related to climatic change by comparing data sets between satellites, as was done in southern Greenland by Zwally and others (Reference Zwally1989) using Seasat and Geosat altimetry.

Acknowledgements

Thanks are extended to R.H. Thomas and W.S. Wilson for supporting this work through the NASA Oceanic Processes Program. Support was also provided by a National Research Council Resident Research Associateship to CS. Lingle at the Oceans and Ice Branch, Code 671, NASA Goddard Space Flight Center. We thank S.N. Stephenson and R.A. Bindschadler for valuable discussions and suggestions, G. Podesta and L.A. Rasmussen for discussions which benefited this work during its early stages, and R.J. Braithwaite and H.H. Thomsen for comments and forwarding results from their field studies in the West Greenland ablation zone. This paper was improved as a result of reviews by D.R. MacAyeal, U.C. Herzfeld, and R.J. Braithwaite.