1. Introduction
Mapping ice-surface elevation with continent-wide coverage is presently best achieved by satellite radar altimetry Since the first satellite that carried an altimeter was already operational in 1978 (Seasat) and has been followed by Geosat (1985–89), the European Remote-sensing Satellites, ERS-1 and ERS-2, and the Japanese Earth Resources Satellite (JERS-1), this type of elevation mapping is useful for monitoring changes in Antarctic and Arctic glaciers, and the average elevation of the Greenland and Antarctic ice sheets. Monitoring changes in ice masses is important for assessing sea-level changes (Reference MeierMeier, 1984) and for investigations and speculations related to global climate change. Accurate three-dimensional models of the actual ice surface serve as a data source for comparison with numerical models and scenarios of ice-sheet breakup or stability (Reference HuybrechtsHuybrechts, 1993). The satellite missions vary in scientific goals, as well as in schedules of repeat times and coverage; geodetic-mission-type data have a dense coverage because the satellite drifts, whereas exact-repeat-type data have a semi-rhombic ground-track pattern with larger gaps between the repeat tracks. All data, however, have gaps that require interpolation prior to map construction. The interpolation is best performed using methods from geostatistics; the theory of regionalized variables (Reference MatheronMatheron, 1963, for general background; Reference Herzfeld, Lingle and LeeHerzfeld and others, 1993, for the method applied to altimetry). There are two approaches for mapping a continent: (a) all-in-one map (e.g. Reference BamberBamber, 1994), or (b) a collection of maps in atlas style. The advantages of method (b) are higher resolution, less distortion due to cartographic projection and higher accuracy. So such maps are more useful for a range of geophysical applications.
There is ongoing work on constructing elevation maps from stereo images of synthetic-aperture-radar (SAR) data, which may yield maps of much higher resolution, but this has yet to be completed. An advantage of using geostatistics for mapping from altimeter data is that a whole set of atlas maps for a given time interval (several months of data are sufficient) can be produced quickly. Time series of maps can be used for monitoring changes of Antarctic ice streams or individual glaciers, or changes in ice-surface elevation averaged over the continent.
An advantage over composite images or datasets from several sources is that a continent-wide map can be derived from a single data source. This is an important aspect in geophysical modeling of ice sheets and glaciers, and for the use of our digital terrain models (DTMs) in other (e.g. climate) models, because the error is homogeneous for the entire area (for error analysis, see Reference Herzfeld, Lingle and LeeHerzfeld and others, 1993).
The maps, or the DTMs, may also be used as base maps and combined with other images or data. An objective of this paper is to demonstrate possible uses of the Geosat and ERS-1 atlas maps. A special focus is the role of the variogram for special-purpose mapping; different techniques are required and best suited for continent-wide atlas mapping, for monitoring changes of the ice surface over a given time span, and for the study of individual glaciers, respectively.
2. Data Acquisition
Altimeter data from Seasat, Geosat, and ERS-1 and 2 are used in mapping the Antarctic ice sheet. Data from an exact repeat mission (ERM) have tracks that repeat exactly at a spacing given by the repeat time, e.g. 35 days during the 1987–89 Geosat-ERM. Data from the Geosat geodetic mission (GM 1995–96) have a high density of ground tracks resulting from the fact that the satellite was allowed to drift.
1995 ERS-1 data are of the geodetic-mission type, the time interval of data used in our maps is 1 February-1 August 1995. Orbital coverage of Geosat altimeter data extends to 72.1° S, ERS-1 coverage to 81.5° S. Data density increases towards the limit of coverage because of the convergence of the ground tracks. Data were processed by the Ice Sheet Altimetry Group at NASA Goddard Space Flight Center (H.J. Zwally J. DiMarzio and coworkers): using Reference Martin, Zwally, Brenner and BindschadlerMartin and others’ (1983) method for retracking; Goddard Earth model (GEM) T2 orbits (Reference MarshMarsh and others, 1989) for data reduction; applying corrections for atmospheric effects and solid Earth tides as described by Reference Zwally, Bindschadler, Brenner, Martin and ThomasZwally and others (1983); slope corrections as described by Reference Brenner, Bindschadler, Thomas and ZwallyBrenner and others (1983); and water-vapor corrections. Our direct link to the Ice Sheet Altimetry Group was employed for transferring ice-data-record (IDR) datasets. From each IDR dataset, those points with retracked and slope-corrected data were retained, prior to mapping. For each map sheet, a track plot is constructed to investigate coverage and ensure that coverage by retracked and slope-corrected data is sufficient. This was the case for all map sheets.
After this processing, "bad" tracks with elevations (a) much lower than the surrounding area, or (b) of about a constant, small (50 m), difference to the surrounding area remained in several ERS-1 maps. An algorithm was developed to identify and remove these bad-track data. Elevation is given with reference to the WGS84 (World Geodetic System) ellipsoid.
3. Atlas Mapping For Continent-Wide Coverage
We use an atlas-mapping scheme, rather than inverting all the Antarctic radar-altimeter (RA) data onto a grid, to produce a single map covering Antarctica (with a hole for the area south of the limit of radar-altimetry coverage). This improves resolution, facilitates mapping of detailed structures, and reduces distortion due to cartographic projection, which is particularly severe for high latitudes and large areas in most algorithms.
An atlas, in the sense of differential analysis (Reference Holmann and RummlerHolmann and Rummler, 1972, p. 63), is a set of maps that (a) covers a given area completely (that is, each point in the area is contained in at least one map); and (b) where projections restricted to areas that appear on two (adjacent) maps (subsets of two maps) are identical at the intersection. In cartography, only property (a) is required for an atlas and the neighborhood relationships need to be matched between sheets.
A useful projection for mapping at high latitude and in atlas form is the Universal Transverse Mercator (UTM) projection (Reference SnyderSnyder, 1987), which results in an orthogonal coordinate system with coordinates in meters. Sufficient overlap of adjacent sheets is convenient for the user of the atlas and necessary to ensure that each point of Antarctica is contained on at least one map despite the distortion of the map edges introduced by the projection algorithm. With the help of a computer program developed especially to facilitate atlas design (Reference Herzfeld, Matassa and MimlerHerzfeld and others, 1999), we decided on the following set-up for projection and map sheets: Map sheets of the ERS atlases are designed to match those of the Geosat atlas. The Antarctic is divided into rows of map sheets (63–68° S, 67–72.1° S, 71°–77° S, 75°–80° S, 78°–81.5° S). In rows 63°–68° S and 67°–72.1° S, maps are of 16° longitude nominal size and overlap 2° on each side, so each map is offset against the next one by 12° of longitude. In rows 71°–77° S, 75°–80° S and 78°–81.5° S, maps are of 36° longitude nominal size, overlap is 6° on each side, offset of two adjacent maps is 24° of longitude. The map names (e.g. m69e61-77n67-721) give the central meridian (69°) and extent of the nominal map area (61°–77° E, 72.1°–67° S (the minus in the north coordinate was dropped in the file names for ease of reading)). The true map area is then the maximal area contained in the nominal map area with the property that it has straight edges in the UTM coordinate system and the grid coordinates are multiples of 3000 (Reference Herzfeld and MatassaHerzfeld and Matassa, 1999). All maps in one row have the same size and grid coordinates.
4. Geostatistical Methods for Elevation Mapping
Interpolation of corrected altimeter-elevation data is carried out using geostatistical methods, in this application ordinary kriging. "Kriging" comprises a family of interpolation/extrapolation methods based on the least-squares optimization principle and usually introduced in a probabilistic framework (Reference MatheronMatheron, 1963; Reference Journel and HuijbregtsJournel and Huijbregts, 1978). However, it can be considered as a numerical-interpolation technique only (Reference HerzfeldHerzfeld, 1992a). Kriging consists of two steps: (1) analysis of spatial structures (variography); and (2) estimation, interpolation and extrapolation (kriging).
4.1. Effect of variogram models on elevation mapping
In the first step, a variogram is calculated according to
where z(xi ), z(xi + h) are measurements at locations xi , xi + h, respectively, inside a region, D, and n is the number of pairs separated by the vector h. The residual variogram is
where
is the drift component. The experimental variograms are calculated in distance classes from the measurements, in our case from the RA data, then are fitted by analytical variogram models. A variogram model describes the type of transition from the strong covariation between closely neighboring samples to weaker covariation of samples farther apart. The variogram model is characterized by its function type, which has to meet the positive-definiteness criterion to ensure existence and uniqueness of the solution of the kriging system. For mapping altimeter data of ice surfaces, a Gaussian variogram model with a relatively high nugget effect is used, because:
-
(1) the ice surface at 3 km resolution is smooth, and the Gaussian model is the model of a mean-square differential process;
-
(2) altimeter data at this resolution have a high signal-to-noise ratio and the spacing of tracks influences the variogram.
The function of the Gaussian variogram model is
The nugget effect, C 0, is the residual variance of resampling in the same location for altimeter data, it also contains contributions caused by surface morphology at a resolution higher than the resolution of the altimeter footprint (called the support). The total sill, C 0 + C 1 is close to the total variance of the data. The concept of a regionalized variable, which is fundamental in geostatistics, perceives that closely neighboring measurements have a higher covariance than measurements spaced farther apart. For distances increasing from the support to the range parameter a in Equation (4), the variogram increases and the regionalization effect decreases; for distances beyond the range, the regionalized variable theoretically behaves like a random variable (probabilistically speaking), and data with a distance larger or equal to the range all have the same influence on the interpolated value (numerically speaking).
For Geosat data
for ERS-1 data
The noisiness of the data is captured by the ratio n 0 of the nugget effect C 0 to the total sill C 0 + C 1 in the following equation:
The nugget-effect ratio, n 0, is 0.422 for Geosat data and 0.705 for ERS-1 data, so ERS-1 data are noisier. The variograms were fitted to data from an area around the grounding zone of Lambert Glacier, a large ice stream in East Antarctica. To avoid edge effects between maps of the same atlas, the same variogram is applied throughout any one atlas.
The effect of using a variogram with a higher nugget-effect ratio is that contour lines are smoother. This is correct, because noisiness in the data does not warrant picking up too many details everywhere.
Selection of the variogram, however, depends not only on the data analysis but also on the mapping purpose. To facilitate comparison, when monitoring changes from one year to another, the same variogram needs to be used for both years. This has been done in the study of Lambert Glacier-Amery Ice Shelf (Reference HerzfeldHerzfeld and others, 1997).
To improve the representation of ice-surface structures in terrain models of individual glaciers, regional variations of the variogram need to be taken into account. To this extent, classes of ice with homogeneous surface properties such as floating ice tongues, grounded glacier ice, steep margins, mountainous terrain, and almost flat inland ice are characterized by specific variogram functions, which may then be used in the interpolation.
4.2. Kriging
The kriging method called "ordinary kriging" (universal kriging of order 0) is better suited for interpolation of RA data than universal kriging of a higher order, because the drift parts modeled by a polynomial component in the higher-order universal kriging methods is likely to create artefacts in the gaps.
In ordinary kriging, the value z 0 = z(x 0) at a node x 0 is estimated by
with data zi = z(xi ) at locations xi (i = 1, …, n) in a neighborhood of the grid node x0. In the Antarctic atlas mapping, the 4 nearest neighbors in each quadrant were used in the interpolation of any given grid node. The coefficients are determined such that the estimation error has minimum variance and the estimation is unbiased, which requires a condition
A solution of the kriging system is obtained using the variogram model specified earlier. (For a derivation of the kriging system, see Reference HerzfeldHerzfeld 1992a, b.) This is carried out for every grid node in the map area. A 3 km grid size is chosen so that the resultant maps or terrain models can be used in glaciodynamic investigations.
5. Examples and Applications
In the sequel, a suite of atlas maps is presented to demonstrate the possible uses and limitations of mapping by geostatistical interpolation from RA data. An analysis of the interpolation is given in Reference Herzfeld, Lingle and LeeHerzfeld and others (1993).
A good example of monitoring changes in an ice stream is presented in the study of Lambert Glacier-Amery Ice Shelf, East Antarctica’s largest ice-stream-ice-shelf system (Reference Herzfeld, Lingle and LeeHerzfeld and others, 1994, Reference Herzfeld1997). Lambert Glacier advanced 10–12 km between 1978 and 1987–89, which corresponds to an advance of about 1 kma−1, as noted by an advance of the grounding zone (Reference Herzfeld, Lingle and LeeHerzfeld and others, 1994).
Applications given here are indexed on the map in Figure 1.
Fig. 1. Map of Antarctica showing location of examples.
5.1. Comparison of radar-altimetry-based map with satellite image: Slessor Glacier
Slessor Glacier is one of the longest coherent ice streams, about 500 km long from Parry Point (the corner point formed by the northern margin of Slessor Glacier and the "coast"-line inland ice/Filchner Ice Shelf) as determined by flowlines visible in satellite imagery (Reference SwithinbankSwithinbank, 1988, p. B100) (Fig. 2). Slessor Glacier and part of the Filchner Ice Shelf are seen on map m333e315-351n78-815 "Filchner Ice Shelf" of the ERS-1 atlas (Fig. 3). The large drainage basin (575 000 km2; Reference SwithinbankSwithinbank, 1988) of Slessor Glacier dominates this map sheet. Slessor Glacier, at its narrowest point (opposite Mount Sheffield in the Shackleton Range), is 30 km wide. It descends west with a constant slope of 0.32° to a flatter area before joining the Filchner Ice Shelf. Surface elevations, 10–30 m higher than the surrounding ice shelf, mark its protrusion into the ice shelf (at 400,000 East,–8,900,000 North).
Fig. 2. Landsat image ofSlessor Glacier (from Reference SwithinbankSwithinbank 1988, p. B100).
Fig. 3. ERS-1 atlas map m333e315-351n78-815 "Filchner Ice Shelf", showing Slessor Glacier and its drainage basin.
In the south, the glacier is bordered by the Shackleton Range. The highest peak of the Herbert Mountains, in the central part of the Shackleton Range, is Mt. Absalom (1642ma.s.L). The Shackleton Range rises to about 1800 m a.s.l. (Reference AlbertsAlberts, 1995). The map kriged from altimeter data reaches over 1600 m in the Herbert Mountains, so the averaging effect that must be expected from both satellite altimetry and kriging is not severe here.
The Landsat image (Fig. 2) shows flowline indicative of rapid motion on the glacier and on the part that is in the Filchner Ice Shelf. Surface roughness in the upper part of Slessor Glacier (1000–1500 ma.s.L; Reference SwithinbankSwithinbank, 1988, p. B99) is characterized by escarpments trending across the flowline, spaced about 10 km apart and attributed to subglacial ridges. Isolated smooth undulation (visible in the Landsat image between Parry Point and Blaiklock Glacier) may be indicative of isolated grounded areas in otherwise floating ice, so this area may be near the grounding line. Comparing this location to the atlas map, the grounding line would coincide with the location of the extension of the line between the inland ice and Filchner Ice Shelf.
5.2. Comparison of Geosat (1985–86) and ERS-1 (1995) atlas maps: Fimbul Ice Shelf
On the one hand, monitoring of ice advance or retreat from an evaluation of RA data, at two points in time and from different satellites, requires knowledge of the general offset, in this case of Geosat and ERS-1. On the other hand, the accuracy of mapping from a given altimeter instrument may be assessed from contoured maps, as demonstrated in the following example. Geosat (1985–86) and ERS-1 (1995) atlas maps m3wellw-5n67-721 "Fimbul Ice Shelf" are given in Figure 4a and b, respectively. The large feature is the ice tongue of Jutulstraumen extending straight north into Fimbul Ice Shelf. Jutulstraumen follows a large geologic trough believed to be of structural origin (Reference Decleir, van Autenboer and CraddockDecleir and Van Autenboer, 1982; after Reference SwithinbankSwithinbank, 1988).
Fig. 4. Maps ofFimbul Ice Shelf (m3wellw-5n67-721). Scale 1:5000000 in atlases, here approximately 1:8700000 due to reduction, (a) Geosat 1985–86 atlas map, contoured from 3 km DTM, calculated from kriging of RA data with a Gaussian variogram model with nugget effect 250 m 2, total sill 593 m 2, nugget-effect ratio 0.422 and range 18 000 m (Geosat variogram). (b) ERS-11995 atlas map, contoured from 3 km DTM, calculated from kriging of RA data with a Gaussian variogram model with nugget effect 43 m2, total sill 61m2, nugget-effect ratio 0.705 and range 16000 m (ERS-1 variogram). (c) Geosat 1985–86 map, contoured from 3 km DTM, calculated from kriging of RA data with ERS-1 variogram. (d) ERS-11995 map, contoured from 3 km DTM, calculated from kriging of RA data with Geosat variogram.
Contours of the Geosat map are much smoother than those of the ERS-1 map, possibly due to more data and higher noise in the latter. Geosat data are prone to noise offshore, caused by signals reflecting off islands and erroneously recorded or corrected. This problem has been largely eliminated from the ERS-1 data. Since the Geosat and the ERS-1 atlases have been calculated with different variograms (cf section 4.1), the ice shelves are not immediately comparable from the atlas maps.
The variogram used for the ERS-1 atlas has a larger smoothing effect than that used for the Geosat atlas, as quantified by the nugget-effect ratio (cf. section 4.1). To study the effects of these two variograms on the cartographic output and to eliminate it from a comparison of different years, the Geosat data have been rekriged using the ERS-1 variogram, and vice versa (cf. Fig. 4c and d). Apparently, the effect of the different variograms is fairly small, the Geosat map still is smoother than the ERS-1 map, and very few features changed after the variogram was altered. Prior to a direct comparison of ice elevation, however, the general offset between Geosat and ERS-1 needs to be known. On all maps, the ice-shelf edge is determined by two near-parallel contour lines (10 m difference) close to each other. The coastline inland ice/Fimbul Ice Shelf has many indentations, probably due to sampling because this error is also found on other ERS-1 maps.
Next to Jutulstraumen, the smaller Schytt Glacier may be identified (490,000 East; –7,950,000 to –7,880,000 North) (cf.Reference SwithinbankSwithinbank,1988,p.B95).
5.3. Detail study of an individual glacier: Williamson Glacier, Sabrina Coast
Williamson Glacier flows into Moscow University Ice Shelf, east of Law Dome, as seen on the ERS-1 atlas map mll7el09-125n63-68 "Sabrina Coast" (Fig. 5a). The details in the DTM may be investigated, as shown in the map of Williamson Glacier in Figure 5b The resolution of the detailed map is the same (3 km grid) as that of the entire atlas. This demonstrates that more information is contained in the DTM than is obvious from a plot of the maps at atlas scale. The models may be used for geophysical study, or for monitoring smaller outlet glaciers of the Antarctic inland ice. Not much information, other than satellite data, is available for this part of East Antarctica.
Fig. 5. (a) ERS-11995 atlas map mll7el09-125n63-68, "Sabrina Coast" (shows location of Williamson Glacier). Scale 1:5 000 000 in atlases, here approximately 1:8 800 000. (b) Detail map of Williamson Glacier. Note that this map is contoured from same 3 km grid as the atlas map.
Acknowledgements
This work was supported in part by the U. S. National Aeronautics and Space Administration Office of Polar Programs under grants NAGW-3790 and NAG5-6114.
