A SAR record of early 21st century change in Greenland

Abstract

Glaciers in Greenland are changing rapidly. To better understand these changes, we have produced a series of seven synthetic aperture radar (SAR) backscatter mosaics for seven winters during the period from 2000 to 2013. Six of the mosaics were created using RADARSAT Fine-Beam data and the seventh used ALOS-PALSAR Fine-Beam Single-Polarization data. The RADARSAT mosaics are radiometrically calibrated and capture changes in the backscatter coefficient related to melt and other events, particularly the strong melting in the summer of 2012. Comparison of features in the ascending-orbit ALOS mosaic and the descending-orbit RADARSAT mosaics indicate that in areas of smooth to moderate topography their locations are consistent to within a few tens of meters. The locations of features identifiable in the RADARSAT mosaics, which were collected with the same imaging parameters, generally agree to within better than the 20 m posting of the data. With such geometric accuracy, these data establish a record of change in Greenland for the early part of the 21st century, thus providing a baseline that can be compared with new radar and optical datasets.

1. INTRODUCTION

Numerous studies have documented substantial change of the Greenland Ice Sheet over the past two decades. Many of the ice sheet's outlet glaciers have sped up (Rignot and Kanagaratnam, 2006; Joughin and others, 2010; Enderlin and others, 2014) with an average speedup for glaciers in the northeast and southwest of >30% (Moon and others, 2012). In particular, the speeds of three of Greenland's largest outlet glaciers, Jakobshavn Isbrae (Joughin and others, 2004, 2014), Helheim glacier (Howat and others, 2005) and Kangerdlugssuaq glacier (Luckman and others, 2006), have varied substantially (Howat and others, 2011). At the same time, surface mass balance (SMB) has declined, largely due to increased summer melt (van den Broeke and others, 2009) and decreased precipitation along the western margin (Seo and others, 2015). As a consequence of increased discharge and reduced SMB, Greenland now contributes ~1 mm to sea rise annually (0.73–1.05 mm a⁻¹) (Shepherd and others, 2012; Enderlin and others, 2014). In addition to speeding up and thinning, many outlet glaciers’ calving fronts have retreated by several hundred meters or more since the early 1990s (Moon and Joughin, 2008; Box and Decker, 2011; Seale and others, 2011; Walsh and others, 2012).

Radar imagery provides an important complement to the growing inventory of high-quality optical imagery collected over the ice sheets (e.g. via the Landsat series of satellites). The usual advantages of synthetic aperture radar (SAR) imaging, such as all-weather, day/night imaging with consistent geometry and illumination, apply to the ice sheets. As with other radars (Nghiem and others, 2012), SAR imagery provides more detailed information about the changes in various ice-sheet snow facies (bare-ice through the dry- snow zone) that are not evident in optical imagery (Fahnestock and others, 1993). High-resolution SAR imagery also provides year-round information about supraglacial lakes, especially whether they remain filled throughout the winter (Poinar and others, 2015).

Through the ongoing Greenland Ice Mapping Project (GIMP), we are producing satellite-derived datasets to document the changes occurring in Greenland (Joughin and others, 2010; Howat and others, 2014). Much of this work aims to measure changes in the flow speeds of the ice sheet's outlet glaciers using interferometric and speckle-tracking techniques applied to SAR data (Joughin and others, 2010, 2012; Moon and others, 2012). Although our major focus has been on velocity, the SAR imagery also provides an important record of overall change in Greenland. Thus, we have produced six ice-sheet wide image mosaics (winters of 2000/01, 2005/06, 2006/07, 2007/08, 2008/09 and 2012/13) using RADARSAT C-band and a single mosaic (winter 2009/10) using ALOS-PALSAR L-band SAR data. The major focus of this paper is to describe the production of these data and to demonstrate their utility.

2. METHODS

The first comprehensive SAR mapping of the Greenland ice sheet was made with the ERS-1 satellite in 1992 (Fahnestock and others, 1993). Five years later, the Canadian Space Agency completely mapped Antarctica with radar for the first time (Jezek, 1999), at which time Greenland was also imaged. These 2000/01 RADARSAT data along with acquisitions from 2005/06 to 2006/07 were used to map change in glacier terminus positions (Moon and Joughin, 2008). These mosaics, however, were uncalibrated and terrain corrected with a relatively coarse DEM. Here we describe the methods used to produce an expanded dataset (seven winters) with greatly improved processing, geometric accuracy and radiometric calibration (RADARSAT only).

A single, complete imaging of Greenland using the fine-beam modes of RADARSAT and ALOS-PALSAR requires image acquisitions along several dozen adjacent, overlapping tracks. To facilitate processing, most of these coast-to-coast swaths were further divided into several-hundred-kilometer-long segments, resulting in more than a hundred images needed to produce a single mosaic of Greenland. Here we describe the steps in the production of the full mosaics.

2.1. SAR data

Our source data were delivered in the form of unfocused (Level-0) SAR image data from the Alaska Satellite Facility (ASF). The RADARSAT image data were all collected using the Fine-Beam-1 mode with HH polarization with a nominal incidence angle range of 37°–41°. All of the ALOS-PALSAR data were collected with Fine-Beam Single-Polarization (HH) mode with a nominal incidence angle of 36.4°–40.8°. We used the Gamma Modular SAR Processor (MSP) (Werner and others, 2000) to process these data into uncalibrated, detected SAR images referenced to a fixed-gain setting with no antenna-pattern correction applied at this stage.

As a spaceborne radar's several-kilometer-wide beam sweeps past a target on the ground, the Doppler frequency of the energy backscattered by the target changes with time (or equivalently with relative position in the beam). To focus the data, a SAR processor correlates (match filters) each point in the output image with its expected Doppler history, which rejects energy from other points, which have different Doppler histories. This process synthesizes an aperture with resolution of a few meters from data collected by a real aperture with resolution of a few kilometers. This procedure requires knowledge of the parameters that describe the expected Doppler history, including the Doppler centroid (center of the returned Doppler spectrum). If the Earth did not rotate, then the Doppler centroid would be zero for a beam that is directed broadside (i.e. orthogonal to the satellite track). While the RADARSAT antenna is nominally directed broadside, its Doppler centroid drifts due to the additional latitude-dependent contribution of the Earth's rotation to the Doppler centroid, complicating processing, which often assumes a constant centroid. Thus, we processed the data with a modified version of the Gamma processor that uses a fixed reference Doppler centroid, but which tracks the drift of the actual Doppler bandwidth about the true centroid (Joughin, 2002). This allows the much longer segments of raw data to be processed into single images. The ALOS-PALSAR data were delivered as individual, overlapping frames, which we concatenated together for processing as longer (several hundred km) segments. Because the PALSAR antenna is nominally steered to zero Doppler, there is little variation of the Doppler centroid along track. As a result, we estimated a Doppler centroid for each scene, which was subsequently used to process the entire image.

2.2. Mosaicking

We produced SAR-image mosaics on a polar-stereographic grid with standard latitude of 70° N and the x-axis extending along 45° E, using the WGS84 as the reference ellipsoid. The MSP SAR processor produces the detected (power) images in deskewed range-Doppler coordinates. Thus, the mosaicker must map between the coordinates of the output grid and the SAR images. To do this, the algorithm first uses the satellite-ephemeris data (state vectors) to determine where a point in the output grid would lie in the radar image if the Earth were perfectly ellipsoidal (i.e. with zero surface elevation everywhere) (Joughin, 1995). Except for the points truly on the ellipsoid, this position is incorrect because the slant-range SAR coordinate depends on the surface elevation relative to the ellipsoid. To correct the positions, the algorithm applies an elevation-dependent range correction to determine the actual radar coordinates for the point in the output grid (Curlander and Mcdonough, 1991). The GIMP-DEM provides the ellipsoidal heights required for this correction (Howat and others, 2014).

In creating a mosaic, as the software cycles through the radar images, it determines the subregion of the output grid that corresponds to the current image. For each of these subregions, the algorithm loops over each point, interpolating a value from the source image where valid data are present. At this stage, the mosaicker uses the GIMP-DEM surface elevation to determine and apply a look-angle-dependent correction for the radar antenna pattern. For RADARSAT, the mosaicker uses a pattern supplied by the ASF and for ALOS-PALSAR it uses the pattern supplied with the Gamma MSP software.

There are many areas where images overlap, so that data from two or more images correspond to the same point in the output grid. Such overlap can produce visible seams in the mosaic due to differences in the images caused by changes in the surface backscatter (energy scattered back toward the radar) between passes (e.g. snowfall or melt events), surface displacement (e.g. fast moving glaciers) or processing artifacts (e.g. imperfect correction for the antenna pattern). To reduce the appearance of such seams, the mosaicker uses a feathering algorithm. First, a distance transform is applied to each image to determine the distance in pixels to the nearest edge, d _e. Next, the algorithm uses these distance values to compute a weighting function that varies linearly from 0 to 1 as a function of d _e, which saturates at 1 for distances greater than a fixed feathering length, d _f. This weighting function is used to scale the images and both the weights and scaled data are summed in output buffers. When all the data have been incorporated into the mosaic, the summed weights are used to scale the result to produce a weighted average as the final result. Points on fast moving glaciers features can move substantially between acquisitions of images in regions of overlap where feathering is applied, producing ‘blurring’ artifacts. While we could have avoided such artifacts by not feathering and using one image or the other in the overlap regions, this would have produced ‘image shear’ artifacts, which also are undesirable. Ultimately, we selected the feathering procedure because it produced better results over slow moving regions, which comprise the majority of area imaged.

2.3. RADARSAT radiometric calibration

The MSP processor produces uncalibrated SAR images. Calibrated SAR images typically are distributed as radar backscatter coefficient σ ⁰, values, which are the ratios of the reflected power per unit area to the power incident on the surface (El-Darymli and others, 2014). Since many SAR processors produce imagery with no knowledge of the topography, σ ₀ typically is computed as though the Earth were flat (i.e. zero elevation relative to the ellipsoid). Because ground-range resolution varies strongly with surface slope, the true σ ₀ has a slope-dependent component. Hence, for a surface with uniform scattering properties (i.e. constant σ ⁰), the radar derived σ ⁰ will appear brighter for slopes facing the radar, while it will appear darker for slopes facing away. Although the mosaicking algorithm uses a DEM, we have chosen not to perform a slope correction. This avoids potential spurious variations in backscatter due to errors in, and limited resolution of, the GIMP DEM, and is consistent with the way the RADARSAT Antarctica Mapping Project and Modified Antarctic Mapping Mission (MAMM) mosaics were produced (Jezek, 1999; Jezek and others, 2003). For the relatively flat interior of the ice sheet, this approach produces values relatively close to the actual values. For the mountainous areas at the coast, it produces an image in which it is easier to distinguish topographic features (i.e. it adds a shaded-relief-like effect).

After application of the antenna pattern correction, the imagery was calibrated as

(1) $${\sigma ^{0}} = a \cdot DN - b.$$

Here the antenna-pattern-corrected output of the processor is given by a digital number (DN), and a and b are processor-dependent calibration coefficients, which were initially unknown. Because σ ⁰ values have a large dynamic range, the backscatter coefficients are provided in dB (i.e. 10 log₁₀[σ ⁰]).

To determine the backscatter coefficients, we calibrated the processor output using the MAMM mosaic of Antarctica (Jezek and others, 2003) to compute the parameters in Eqn (1). The MAMM mosaic was calibrated using processor parameters determined from global surveys of targets with known backscatter properties, and is thus a good reference surface to calibrate our processor. To do this, we processed uncalibrated RADARSAT images for 17 areas around Antarctica, which we multi-looked (averaged) 4 pixels in range and 6 pixels in azimuth. After correction for the antenna pattern, we extracted data from areas of overlap with the MAMM mosaic, smoothing both datasets to 400 m resolution. For each area we performed a least-squares fit to extract the calibration parameters. Based on this analysis, we used the mean values of a = 0.03663 and b = 0.0058 to calibrate the Greenland RADARSAT mosaics. With these values the mean difference between the MSP processed data and the MAMM mosaic was 0.0 with a standard deviation of 0.7 dB. Some of this variance likely reflects the fact that we could not ensure that we used the same images as those that went into the MAMM mosaic.

We used a single set of parameters and antenna pattern to calibrate the RADARSAT data. As a consequence, potential variability in the instrument performance or antenna pattern with time could affect our radiometric precision (Srivastava and others, 2007). To assess the extent to which such drift might be present in our data, we examined average radar backscatter coefficients for three 12 km × 12 km regions of exposed rock where we expect natural variation in backscatter to be small. For these regions, the data were consistent to within ~1.5 dB with a decreasing trend of ~0.11 dB a⁻¹ (−0.09, −0.13, −0.12 dB a⁻¹ at three sites). We did not attempt to correct for this apparent drift. The MAMM mosaic that we used as our reference was collected at the same time as the 2000/01 Greenland mosaic. Hence, the 2000/01 mosaic should have the least uncertainty in an absolute sense in so much as the reference MAMM mosaic is correct.

Because we produced only a single mosaic and did not have a calibration reference for ALOS-PALSAR, we did not calibrate the L-band data. Nonetheless, as an image product it is still useful for studying changes such as glacier retreat (see examples below).

2.4. Geometric accuracy

Geolocation errors of 10s to 100s of meters are often introduced by errors in the satellite-along-track time relative to the state-vectors, δ _a, and errors in the range-delay times, δ _r. Thus, the SAR images also require geometric calibration to improve location accuracy.

In our early processing of the RADARSAT mosaics, we noticed along-track location inconsistencies of up to ~70 m between the adjacent tracks. To fix this, we cross-correlated patches in adjacent, overlapping images of the 2005/06 data at 15 m resolution to determine the relative timing errors. We used these relative values of δ _a to adjust the satellite timing to bring all the tracks into relative agreement and then corrected by an additional term so that the mean shift (for 54 data takes) was zero. For the other C-band mosaics, we used a similar registration procedure except that we determined values of δ _a by aligning the data for each track with the 2005/06 mosaic. The ALOS-PALSAR along-track positions were consistent from track to track, but required a uniform correction of δ _a = 0.062 s to align the mosaicked results with independent data (e.g. LandSat and TerraSAR-X).

Comparison of the RADARSAT data with other independent images indicated that there was an approximately constant (for a given year) range delay. To fix this error, we located several easily identifiable features that were visible in both the 2005/06 mosaic and other well-geolocated data (e.g. WorldView and TerraSAR-X). Using these features, we determined that a correction in δ _r equivalent to 70 m was required to align the 2005/06 RADARSAT data with independent data. After evaluating all the RADARSAT data, we used a fixed δ _r correction for each year ranging from 66 to 70 m. For the ALOS data no range correction was required.

After application of the range and along-track timing corrections, all the RADARSAT mosaics were internally consistent to within better than the 20 m pixel spacing of the final mosaics. Comparison with other data (Landsat, WorldView, TerraSAR-X) indicates that mean errors are also <20 m. Thus, the remaining source of error is the GIMP DEM used for terrain correction (Howat and others, 2014). For a nominal incidence angle of 38°, the horizontal range-dependent location error at a point should be a factor ~1.25 (1/tan[38°]) greater than the corresponding elevation error (Curlander and Mcdonough, 1991). For the ice sheet the GIMP-elevation error is 8.5 m (Howat and others, 2014), indicating subpixel (<20 m) location errors for ice-covered regions. For ice-free terrain the elevation errors assessed in different areas of Greenland vary substantially (8–40 m), indicating location errors of up to 50 m. In areas with severe slopes such as mountainous regions along the southeast coast (see Fig. 6 in Howat and others, 2013), the DEM errors may be substantially larger.

In general, the DEM-dependent errors are common to all RADARSAT mosaics, so relative positions are consistent to better than 20 m. An exception occurs where glaciers have thinned by tens of meters over the 12 a period during which the images were acquired. In such regions, horizontal location errors are similar (~1.25×) in magnitude to the corresponding elevation change.

While the RADARSAT mosaics are relatively self-consistent, location differences between the RADARSAT mosaics and ALOS-PALSAR mosaic can be substantially larger. Such differences arise because the RADARSAT images were acquired along descending orbits (satellite moving toward the equator in the northern hemisphere) and the ALOS-PALSAR data were acquired along the ascending orbits (satellite moving toward the pole). As a result, points in the mosaics were imaged from opposing sides by the two satellites. Thus, if the location error due to terrain distortion is ΔX in a RADARSAT mosaic, the corresponding error in the ALOS-PALSAR image will be roughly in the opposite direction in the ALOS-PALSAR image, so the relative error between the images will be $\sim 2{\rm \Delta} X$ , assuming that the elevation error is the same for both images. In addition, because of this opposite-side viewing, areas that are shadowed in the RADARSAT image will be laid-over or foreshortened in the ALOS-PALSAR image and vice versa. Foreshortening occurs when radar-facing slopes are such that slant-range pixels project to a large area on the ground, causing the radar-side of a mountainous feature to be compressed in slant-range coordinates. Layover occurs when slopes are such that points that are farther in ground range are nearer when mapped into slant-range coordinates than other points that are actually closer in ground range (i.e. a peak may be closer to the radar due to its height, even though it is farther way horizontally than a lower elevation point). Although the terrain correction can fix much of this distortion, differences can be large (tens of meters) in areas of extreme slopes.

2.5. Resolution

We produced all the mosaics with postings of both 20 and 100 m. All the RADARSAT data were collected using the FN1 (Fine-1) beam, which has a single-look resolution of ~4.6 m in the slant-range direction and ~5.6 m in the azimuth direction. For the 20 m mosaics, the data were multi-looked by incoherently averaging by 3 pixels in range by 4 pixels in azimuth to reduce speckle. This degree of averaging translates into a nominal 22 m×22 m ground resolution. Actual resolution varies with incidence angle (Curlander and Mcdonough, 1991). For the 100 m mosaics, we multi-looked (downsampled by averaging groups of adjacent n×m pixels) the data using 12 pixels in range and 18 pixels in azimuth, which yields a nominal 89 m×99 m ground resolution. The 20 m ALOS-PALSAR mosaics were produced from images multi-looked 3 pixels in range by 6 pixels in azimuth, which yields a nominal ground resolution of 22 m×19 m. Finally, the 100 m L-band mosaics were multi-looked 12 pixels in range by 24 pixels in azimuth (~88 m×86 m).

3. RESULTS

The main purpose of this section is to illustrate the quality and properties of the data. Thus, we have selected several examples to illustrate both the quality and the potential applications of the data. We have not performed a comprehensive analysis of the data since our main goal is to distribute a well-documented dataset to the broader community for further analysis. All mosaics are available from the National Snow and Ice Data Center (NSIDC) (https://nsidc.org/data/nsidc-0633).

3.1. Complete mosaics

Because no set of data acquisitions from a single year provided complete SAR coverage of Greenland, we also produced an as-complete-as-possible SAR mosaic at each frequency using data collected over multiple winters. The resulting mosaics are shown in Figure 1. The C-band mosaic (left panel Figure 1) is largely based on data collected over the 2005/06 winter. We filled gaps in southern Greenland using data from other years. One swath from 2009 (orbit 68869) was substantially brighter than nearby swaths (possibly an error in the reported automatic gain control values or other parameter in the metadata), so we scaled it by a factor of 0.72 to produce a better match to the surrounding images. The ALOS-PALSAR mosaic was created almost entirely from data collected in the 2009/10 winter. The few L-band images from other years that were used to fill gaps were mostly over ice-free coastal areas.

Fig. 1. Mosaics of SAR data collected with (left) RADARSAT operating at C-band and (right) ALOS-PALSAR operating at L-band. Images are uncalibrated and have been nonlinearly stretched to enhance contrast. Each mosaic contains data from multiple years to provide nearly complete coverage. White, red and blue boxes indicate locations of results shown in Figures 3 to 5, respectively.

3.2. Calibrated annual mosaics

Figure 2 shows the calibrated RADARSAT mosaics for each winter during which data were acquired. We use the term winter loosely here to describe the period when there is little or no surface melting. Thus, a few images were acquired as early as late September and a few as late as April. Nonetheless, for the entire set of RADARSAT mosaics only ~1% of the data were collected in September and April and ~91% of data were acquired in December through February (a shape file is included with the distributed products to provide the date and orbit information for the source images). Each winter, imaging was repeated up to four times along most tracks. In producing the mosaics, we picked a single-coverage (neglecting overlap) subset of these images that provided the maximum coverage, while attempting to minimize the period over which the data were collected (i.e. the images are clustered in time as closely as possible). Most of the gaps in coverage are in the south, where the instrument collected data in other modes for sea-ice mapping in shipping lanes.

Fig. 2. Calibrated RADARSAT σ ⁰ mosaics collected over six winters in the period from 2000 to 2013. Images were collected over the intervals: September 21, 2000 to January 23, 2001; December 24, 2005 to April 4, 2006; December 30, 2006 to February 4, 2007; November 22, 2007 to March 30, 2008; January 10, 2009 to February 5 2009; January 15, 2013 to March 24, 2013.

As detailed by Fahnestock and others (1993), backscatter levels correspond roughly to the various snow facies on the ice sheet (Cuffey and Paterson, 2010). This correspondence is largely determined by how seasonal melt influences the scattering properties. In the central part of the ice sheet, particularly before 2012/13, the dry-snow zone is relatively dark (~−16 to −19 dB) due to weak volume scattering from fine-grained snow and layering that has not been modified by melt. In the percolation-zone, refreezing of seasonal melt produces ice lenses, pipes and layers, which backscatter strongly (~ −3 to −4 dB). Although there is greater melting in the wet-snow zone, the refreezing is more uniform, producing less density variation and backscatter (~−5 to −6 dB) relative to the percolation zone. Finally, recent snowfall is nearly transparent, so in the bare-ice zone most of the radar return is caused probably by relatively weak backscatter (~−13 to −15 dB) from the rough ice surface exposed the prior summer. Additional brightness variations are introduced by scattering from features such as stream channels, frozen-over lakes and crevasses.

While Figure 2 shows the full time series, it is difficult to see interannual variability in backscatter. To better demonstrate this variability, Figure 3 shows the difference between σ ⁰ for each year and the mean backscatter coefficient, $\overline {{\sigma ^0}} $ , for the relatively stable period from 2000 to 2009 (we excluded 2012/13 from the average since it is such an outlier). The results indicate minor (a few dB) variations in backscatter through 2009. The radar backscatter from the northeast region of the ice sheet is low (~2–4 dB below average) in 2000/01. This region brightened as of 2005/06 and the backscatter remained relatively stable for the next 2 a. By the winter of 2008/09, much of the ice-sheet margin brightened by ~1–3 dB. The high-elevation ice-sheet interior brightened substantially (~5–10 dB), following the extreme melt event in summer 2012 (Nghiem and others, 2012). This brightening is due to the refreezing of summer melt to form ice layers, lenses, layers and pipes similar to those in the percolation zone (i.e. at least transiently, the percolation zone effectively extended over the entire interior of the ice sheet in 2012).

Fig. 3. Difference RADARSAT σ ⁰ values relative to the average for 2000 through 2009, $\overline {{\sigma ^{0}}} $ . Results are only plotted for areas where there are data for all the years.

To examine interannual changes in backscatter in more detail, the left-hand panels in Figure 4 show a region in Southwest Greenland that extends, west to east, from the ice-sheet margin to the lower percolation zone and includes the region covered by the K-transect (van de Wal and others, 2008). The bright region around the margin of the ice sheet approximately indicates the extent of the percolation zone (Fahnestock and others, 1993). The most notable feature is the 2012/13 inland shift by several kilometers of the brightest region, indicating the upslope migration of the percolation zone during the extreme melt event. Also evident in the sequence of images is the migration of supraglacial lakes to higher elevations on the ice sheet (Howat and others, 2013).

Fig. 4. (left) Calibrated C-band-image examples showing an area along the western margin of the ice sheet (see white box in Fig. 1 for location). White (2000/01) and black (2012/13) lines in top and bottom panels indicate the gradation in radar brightness that roughly indicates the transition from wet snow to the percolation zone. (right) Detail of the area around a supraglacial lake indicated by white box in the 2000/01 image.

To better illustrate the detail captured by the full-resolution 20 m images, the right-hand panels of Figure 4 show the area around a supraglacial lake at an elevation of 1325 m (WGS84). At this resolution features such as large outflow streams from the lake are visible. Here, the bright surface on the lake during some years (2000/01, 2005/06 and 2007/08) likely indicates a thick ice cover over water remaining from a lake that did not drain completely the prior summer (Jeffries and others, 1994). In other years (2006/07 and 2008/09), there appears to be an ice cover, but the backscatter is substantially darker. This could indicate wet snow on the surface or perhaps a lake that drained after the ice cover formed. In 2012/13, there is no sign of ice cover and a feature likely corresponding to a moulin is visible, together indicating this lake drained during the prior summer.

3.3. Change in ice-sheet geometry

As noted above, except where there are large changes in surface elevation, the C-band image mosaics are co-registered to within better than the 20 m posting. As an illustration of this level of precision, Figure 5 shows an example of an unnamed, small, shrinking ice cap in Northern Greenland (see Figure 1 for location). We digitized the outlines of the ice cap in the 2000/01 and 2012/13 images and plotted them over the images from all years. These results show contraction of the ice-cap margin over the 12 a period, varying along the margin over a range from ~100 to 500 m, which is well resolved with the 20 m posting. The retreat is relatively steady, as indicated by the margin positions in the intermediate years. Because the ice cap is thinning, geolocation accuracy is affected by our use of a single DEM for all years. The slopes of these ice caps are relatively low (~0.025–0.05); however, yielding errors of ~3–6%, which is comparable with the digitization error (~1 pixel). Close examination of the mosaics indicates that many of the other small ice caps in Northern Greenland are retreating similarly.

Fig. 5. An example showing shrinkage of a small ice cap in Northern Greenland (see red box in Fig. 1 for location). The lines show the hand-digitized margins for 2000/01 (red) and 2012/13 (blue).

The radar mosaics provide a nearly complete sampling for the winter terminus positions of Greenland's outlet glaciers (Moon and Joughin, 2008). Figure 6 provides an example showing a time series of images for three glaciers along the west coast of Greenland (see Fig. 1 for location). In this example, the northernmost glacier retreated by ~1 km from 2000/01 to 2005/06, where it maintained an approximately stable position for at least the next 4 a. Over the period from 2008/09 to 2012/13 the glacier retreated ~3 km farther inland. By contrast, the calving fronts of the two southern glaciers remained steady with interannual variability of <300 m.

Fig. 6. Examples of glacier terminus retreat and stability over the 12 a period spanned by the mosaics (see blue box in Fig. 1 for location). Images are uncalibrated with a stretch to enhance contrast. Color lines show the terminus position in 2000/01 (red), 2005/06 (green), 2006/07 (orange), 2007/08 (yellow) and 2008/09 (magenta). From top to bottom the glacier names are Sermeq Silarleq, Kangilleq and Sermilik.

The examples in Figures 5 and 6 reveal changes mapped using images acquired with the same imaging geometry, along the descending portions of the RADARSAT orbits. The 2009/10 ALOS mosaic, however, was created using the data acquired entirely along ascending tracks. As noted earlier, evaluation of ascending and descending images provides a worst-case comparison since the surface is imaged from opposite sides, causing terrain-induced geolocation errors in opposing directions. To examine the types of differences that arise, Figure 7 shows a comparison of the 2008/09 C-band mosaic (descending tracks) with the 2009/10 L-band mosaic (ascending tracks). In areas with steep slopes, there is strong layover or shadowing (Figure 7 top), introducing distortion that cannot be fully corrected even with a relatively high-resolution DEM. On flatter, stationary areas, the locations of features tend to line up well, generally to within ~50 m or less.

Fig. 7. Comparison of (left) C-band from 2008/09 and (right) L-band imagery from 2009/10. (top) Same areas as Figure 6 with ice-front locations as indicated in that figure with the addition of the 2012/13 front in blue. (bottom) Small ice-cap example from Figure 5. Color lines show the terminus position/ice-cap margin in 2000/01 (red), 2005/06 (green), 2006/07 (orange), 2007/08 (yellow), 2008/09 (magenta), 2009/10 (cyan) and 2012/13 (blue).

For the two southern glaciers in Figure 7, the L-band terminus positions fall well within the range of the interannual variability indicated by the C-band data. For the northernmost glacier, the L-band image indicates a little (<1 km) of retreat from 2008/09 to 2009/10, indicating that the 3 km retreat from 2008/09 to 2012/13 (Figure 6) likely did not occur in a single year. Thus, despite the differences in imaging geometry, in this example comparison of ascending and descending images is sufficient to detect changes in terminus position with errors smaller than the natural fluctuations that occur over a typical winter. Errors would be substantially larger for areas with extremely rapid elevation change, such as near the terminus of Jakobshavn Isbrae, where the surface elevation has lowered by >200 m since the 1990s. In the case of Jakobshavn Isbrae, if the DEM was from the 1990s, then the absolute error would be ~250 m for present day positions. When comparing ascending and descending images, this error would be in opposing directions, so that the relative error would be nearly half a kilometer (the GIMP DEM corresponds roughly to the middle of our period, so the differences are smaller).

The bottom of Figure 7 shows a comparison of the small ice cap from Figure 5 for C-band (2008/09) and L-band (2009/10) mosaics both overlaid with the digitized outlines from 2008/09 and 2009/10. The disagreement between the digitized outlines ranges from a few meters up to ~100 m, reflecting a combination of interannual change and digitization and geolocation errors, which are difficult to evaluate individually. Comparison of digitized outlines of the 2007/08 and 2008/09 mosaics with like-imaging geometries (not shown) provides a level of agreement that is similar to that of the ascending/descending mosaics. This agreement indicates that both C- and L-band imagery are well registered in regions with relatively flat terrain and errors in the DEM used for terrain correction are small (a few tens of meters).

4. DISCUSSION

Moon and Joughin (2008) used earlier versions of some of the mosaics shown in Figure 2 (2001/02 to 2006/07) to estimate the retreat of outlet-glacier termini around Greenland. These earlier mosaics were terrain-corrected and geolocated with a 1 km posting DEM, which was derived primarily from radar altimetry (Bamber and others, 2001). While this DEM provides good accuracy in the flat interior, height errors may range from 10s to 100s of meters in topographically rough coastal areas, leading to large (10s to 100s of meters) location errors in the SAR mosaics. These errors are common, however, to images with the same viewing geometry and processed with the same DEM. Thus, earlier estimates of terminus change were relatively unaffected since the largely common component of the geolocation error cancels when differencing positions.

Lack of high-resolution, accurate DEMs has long been an issue in processing well-located SAR imagery, particularly for high-latitude regions not covered by the Shuttle Radar Topography Mission (SRTM) (Farr and others, 2007). Our new mosaics (Figure 2) use the 30 m GIMP DEM (Howat and others, 2014) for terrain correction and geolocation, greatly improving absolute accuracy, especially in regions with steep topography. As a result, even images acquired using substantially different viewing geometries (e.g. ascending and descending) compare relatively well (a few tens of meters difference). Only in regions of heavy layover or shadow (e.g. steep fjord walls) it is difficult to compare images with different geometries. Thus, the L- and C-band image mosaics provide a good baseline observation of the geometry of the Greenland ice sheet during the first part of the 21st century. If such images could have been produced throughout the 20th century, then we would have a much better context with which to interpret the current rates of retreat.

Although the SAR data span a 12 a period (late 2000 to early 2013), we terrain-corrected the entire dataset using the GIMP DEM, the elevations of which roughly correspond to the middle of this period (Howat and others, 2014). Any temporal inconsistency between the data and DEM produces negligible errors for the exposed bedrock around the margin, which is fixed (neglecting mm a⁻¹-scale uplift). Similarly, the elevation of most of the interior ice sheet is evolving slowly (tens of cm a⁻¹ or less), yielding only small location errors due to changes in surface elevation. A small percentage of the ice-sheet area, however, is thinning rapidly (>5 m a⁻¹ or more). For these regions, our errors scale with the elevation change relative to the reference GIMP DEM. Over the period our data were collected, such changes are generally within the degree of intraannual variability (e.g. <100 m a⁻¹). Ideally, though, DEMs should be regularly updated in areas of rapid change to maintain a high degree of geometric accuracy.

Although the images shown in Figure 1 are not radiometrically calibrated, qualitatively they illustrate some notable similarities and differences with C- and L-band imaging. The dry-snow zone is dark in both mosaics, but generally is darker at L-band, likely indicating reduced volume scattering from snow grains at the longer (23.6 cm) L-band wavelength relative to the shorter (5.6 cm) C-band wavelength. Similarly, in the bare-ice zone, the surface is noticeably darker at L-band, likely indicating less surface backscatter at the longer wavelength. In addition, the penetration depth is greater at L-band (~3 m vs ~1 m) in the bare-ice zone (Rignot and others, 2001). As a consequence, buried reflectors such as crevasses are far more visible at L-band, especially against the generally darker bare-ice background. Because of the greater L-band penetration in more homogeneous ice, the wet-snow zone is substantially darker at L-band. Finally, in the percolation zone the scale of buried melt features is such that they produce a bright return at both C- and L-bands. Thus, although there is some commonality, both frequencies provide unique information about various surface and near-surface processes.

5. CONCLUSIONS

We have presented a new set of SAR mosaics that record change in Greenland from 2000 to 2013. The examples we provide serve to illustrate and document the accuracy and utility of the data rather than provide a complete analysis of the data. The data are freely available online at the NSIDC with the goal that the broader community will explore the dataset to further understand the changes that have occurred in Greenland since the beginning of the 21st century. While the sensors that collected these data no longer exist, there are several new SARs regularly collecting data over Greenland. Hence, the record we have begun can be carried forward to the future.