Journal cover Journal topic
The Cryosphere An interactive open-access journal of the European Geosciences Union
Journal topic
TC | Articles | Volume 13, issue 2
The Cryosphere, 13, 427–449, 2019
© Author(s) 2019. This work is distributed under
the Creative Commons Attribution 4.0 License.
The Cryosphere, 13, 427–449, 2019
© Author(s) 2019. This work is distributed under
the Creative Commons Attribution 4.0 License.

Research article 05 Feb 2019

Research article | 05 Feb 2019

Four decades of Antarctic surface elevation changes from multi-mission satellite altimetry

Four decades of Antarctic surface elevation changes from multi-mission satellite altimetry
Ludwig Schröder1,2, Martin Horwath1, Reinhard Dietrich1, Veit Helm2, Michiel R. van den Broeke3, and Stefan R. M. Ligtenberg3 Ludwig Schröder et al.
  • 1Technische Universität Dresden, Institut für Planetare Geodäsie, Dresden, Germany
  • 2Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany
  • 3Institute for Marine and Atmospheric Research Utrecht, Utrecht University, Utrecht, the Netherlands

Correspondence: Ludwig Schröder (

Back to toptop

We developed a multi-mission satellite altimetry analysis over the Antarctic Ice Sheet which comprises Seasat, Geosat, ERS-1, ERS-2, Envisat, ICESat and CryoSat-2. After a consistent reprocessing and a stepwise calibration of the inter-mission offsets, we obtained monthly grids of multi-mission surface elevation change (SEC) with respect to the reference epoch 09/2010 (in the format of month/year) from 1978 to 2017. A validation with independent elevation changes from in situ and airborne observations as well as a comparison with a firn model proves that the different missions and observation modes have been successfully combined to a seamless multi-mission time series. For coastal East Antarctica, even Seasat and Geosat provide reliable information and, hence, allow for the analysis of four decades of elevation changes. The spatial and temporal resolution of our result allows for the identification of when and where significant changes in elevation occurred. These time series add detailed information to the evolution of surface elevation in such key regions as Pine Island Glacier, Totten Glacier, Dronning Maud Land or Lake Vostok. After applying a density mask, we calculated time series of mass changes and found that the Antarctic Ice Sheet north of 81.5 S was losing mass at an average rate of -85±16 Gt yr−1 between 1992 and 2017, which accelerated to -137±25 Gt yr−1 after 2010.

1 Introduction
Back to toptop

Satellite altimetry is fundamental for detecting and understanding changes in the Antarctic Ice Sheet (AIS; Rémy and Parouty2009; Shepherd et al.2018). Since 1992, altimeter missions have revealed dynamic thinning of several outlet glaciers of the West Antarctica Ice Sheet (WAIS) and have put narrow limits on elevation changes in most parts of East Antarctica. Rates of surface elevation change are not constant in time (Shepherd et al.2012). Ice flow acceleration has caused dynamic thinning to accelerate (Mouginot et al.2014; Hogg et al.2017). Variations in surface mass balance (SMB) and firn compaction rate also cause interannual variations of surface elevation (Horwath et al.2012; Shepherd et al.2012; Lenaerts et al.2013). Consequently, different rates of change have been reported from altimeter missions that cover different time intervals. For example, ERS-1 and ERS-2 data over the interval 1992–2003 revealed negative elevation rates in eastern Dronning Maud Land and Enderby Land (25–60 E) and positive rates in Princess Elizabeth Land (70–100 E) (Wingham et al.2006b), while Envisat data over the interval 2003–2010 revealed the opposite pattern (Flament and Rémy2012). Two large snowfall events in 2009 and 2011 induced stepwise elevation changes in Dronning Maud Land (Lenaerts et al.2013; Shepherd et al.2012).

As a consequence, mean linear rates derived from a single mission have limited significance in characterizing the long-term evolution of the ice sheet (Wouters et al.2013). Data from different altimeter missions need to be linked over a time span that is as long as possible in order to better distinguish and understand the long-term evolution and the natural variability of ice sheet volume and mass.

Missions with similar sensor characteristics have been combined, e.g., by Wingham et al. (2006b, ERS-1 and ERS-2) and Li and Davis (2008, ERS-2 and Envisat). Fricker and Padman (2012) use Seasat, ERS-1, ERS-2 and Envisat to determine elevation changes of Antarctic ice shelves. They apply constant biases, determined over open ocean, to cross-calibrate the missions. In contrast to ocean-based calibration, Zwally et al. (2005) found significant differences for the biases over ice sheets with a distinct spatial pattern (see also Frappart et al.2016). Khvorostovsky (2012) showed that the correction of inter-mission offsets over an ice sheet is not trivial. Paolo et al. (2016) cross-calibrated ERS-1, ERS-2 and Envisat on each grid cell using overlapping epochs, and Adusumilli et al. (2018) extended these time series by including CryoSat-2 data. We use a very similar approach for conventional radar altimeter measurements with overlapping mission periods. Moreover, we also include measurements of the nonoverlapping missions Seasat and Geosat and measurements with different sensor characteristics, such as ICESat laser altimetry or CryoSat-2 interferometric synthetic aperture radar (SARIn) mode, making the combination of the observations even more challenging.

Here we present an approach for combining seven different satellite altimetry missions over the AIS. Using refined waveform retracking and slope correction of the radar altimetry (RA) data, we ensure the consistency of the surface elevation measurements and improve their precision by up to 50 %. In the following stepwise procedure, we first process the measurements from all missions jointly using the repeat-altimetry method. We then form monthly time series for each individual mission data set. Finally, we merge all time series from both radar and laser altimetry. For this last step, we employ different approaches of inter-mission offset estimation, depending on the temporal overlap or nonoverlap of the missions and on the similarity or dissimilarity of their altimeter sensors.

Figure 1Spatial and temporal coverage of the satellite altimetry data used in this study. The colors denote the maximum southern extent of the measurements (dark blue: 72 S, light blue: 81.5 S, orange: 86 S, red: 88 S) and thus the size of the respective polar gap.


We arrive at consistent and seamless time series of gridded surface elevation differences with respect to a reference epoch (09/2010; in the format of month/year) which we made publicly available at (Schröder et al.2019). The resulting monthly grids with a 10 km spatial resolution were obtained by smoothing with a moving window over 3 months and a spatial Gaussian weighting with 2σ=20km. We evaluate our results and their estimated uncertainties by a comparison with independent in situ and airborne data sets, satellite gravimetry estimates and regional climate model outputs. We illustrate that these time series of surface elevation change (SEC) allow for the study of geometry changes and derived mass changes of the AIS in unprecedented detail. The recent elevation changes of Pine Island Glacier in West Antarctica, Totten Glacier in East Antarctica and Shirase Glacier of Dronning Maud Land in East Antarctica are put in context with the extended time series from satellite altimetry. Finally, we calculate ice sheet mass balances from these data for the covered regions. A comparison with independent data indicates a high consistency of the different data sets but also reveals remaining discrepancies.

2.1 Altimetry data

We use the ice sheet surface elevation observations from seven satellite altimetry missions: Seasat, Geosat, ERS-1, ERS-2, Envisat, ICESat and CryoSat-2. Figure 1 gives an overview of their temporal and spatial coverage. The data of the two early missions, Seasat and Geosat, were obtained from the Radar Ice Altimetry project at Goddard Space Flight Center (GSFC). For the ESA missions, we used the data of the REAPER reprocessing project (Brockley et al.2017) of ERS-1 and ERS-2, the RA-2 Sensor and Geophysical Data Record (SGDR) of Envisat in version 2.1 and Baseline C Level 2I data of CryoSat-2. For ICESat we used GLA12 of release 633 from the National Snow and Ice Data Center (NSIDC). Further details concerning the data set versions used and the data-editing criteria applied to remove corrupted measurements in a preprocessing step are given in the Supplement.

As illustrated in Fig. 1, due to the inclination of 108, Seasat and Geosat measurements only cover the coastal regions of the East Antarctic Ice Sheet (EAIS) and the northern tip of the Antarctic Peninsula Ice Sheet (APIS) north of 72 S, which is about 25 % of the total ice sheet area. With the launch of ERS-1, the polar gap was reduced to areas south of 81.5 S, resulting in coverage of 79 % of the area. The polar gap is even smaller for ICESat (86 S) and CryoSat-2 (88 S), leading to a nearly complete coverage of the AIS in recent epochs.

ERS-1 and ERS-2 measurements were performed in two different modes, distinguished by the width of the tracking time window and the corresponding temporal resolution of the recorded waveform. The ice mode is coarser than the ocean mode, in order to increase the chance of capturing the radar return from rough topographic surfaces (Scott et al.1994). While the ice mode was employed for the majority of measurements, a significant number of observations have also been performed in the ocean mode over Antarctica (22 % for ERS-1, 2 % for ERS-2). We use the data from both modes, as the ocean mode provides a higher precision while the ice mode is more reliable in steep terrain (see Figs. S1 and S3 in the Supplement). However, as there is a regionally varying bias between the modes, we treat them as two separate data sets, similar to Paolo et al. (2016).

2.2 Reprocessing of radar altimetry

Compared to measurements over the global oceans, pulse-limited radar altimetry (PLRA) over ice sheets requires a specific processing to account for the effects of topography and the dielectric properties of the surface (Bamber1994). To ensure consistency in the analysis of PLRA measurements, processed and provided by different institutions, we applied our own method for retracking and slope correction.

The slope correction is applied to account for the effect of topography within the beam-limited footprint (Brenner et al.1983). Different approaches exist to apply a correction (Bamber1994), but this effect is still a main source of error in RA over ice sheets. The “direct method” uses the surface slope within the beam-limited footprint to obtain a correction for the measurement in the nadir direction. In contrast, the “relocation method” relates the measurement towards the likely true position up slope. While the direct method has the advantage that the measurement location is unchanged, which allows an easier calculation of profile crossovers or repeat-track parameters, the relocation method has lower intrinsic error (Bamber1994). A validation using crossovers with kinematic GNSS-profiles (Schröder et al.2017) showed that, especially in coastal regions, the direct method leads to significantly larger offsets and standard deviations, compared to the relocation method. Roemer et al. (2007) developed a refined version of the relocation method, using the full information of a digital elevation model (DEM) to locate the point of closest approach (POCA) within the approximately 20 km beam-limited footprint. We applied this method in our reprocessing chain using the DEM of Helm et al. (2014). The CryoSat-2 measurements, used for this DEM, have very dense coverage, and hence very little interpolation is necessary. Compared to the DEM of Bamber et al. (2009), this significantly improves the spatial consistency. We optimized the approach of Roemer et al. (2007) with respect to computational efficiency for application over the entire ice sheet. Instead of searching the POCA with the help of a moving window of 2 km (which represents the pulse-limited footprint) in the DEM-to-satellite grid, we applied a Gaussian filter with σ=1km to the DEM itself to resemble the coverage of a pulse-limited footprint. Hence, instead of the closest window average, we can simply search for the closest cell in the smoothed grid, which we use as the coarse POCA location. In order to achieve a subgrid POCA location, we fit a biquadratic function to the satellite-to-surface distance within a 3×3 grid cell environment around the coarse POCA grid cell and determine the POCA according to this function.

The retracking of the return signal waveform is another important component in the processing of RA data over ice sheets (Bamber1994). Functional fit approaches (e.g., Martin et al.1983; Davis1992; Legrésy et al.2005; Wingham et al.2006b) are well established and allow the interpretation of the obtained waveform shape parameters with respect to surface and subsurface characteristics (e.g., Lacroix et al.2008; Nilsson et al.2015). However, the alternative approach of threshold retrackers has proven to be more precise in terms of repeatability (Davis1997; Nilsson et al.2016; Schröder et al.2017). A very robust variant is called ICE-1, using the “offset center of gravity” (OCOG) amplitude (Wingham et al.1986). Compared to the waveform maximum, the OCOG amplitude is significantly less affected by noise (Bamber1994). Davis (1997) compared different retrackers and showed that a threshold-based retracker produces remarkably higher-precision results (especially with a low threshold as 10 %), compared to functional-fit-based results. We implemented three threshold levels (10 %, 20 % and 50 %) for the OCOG amplitude, which allowed us to analyze the influence of the choice of this level, similar to Davis (1997).

In addition to PLRA, we also use the SARIn mode data of CryoSat-2, reprocessed by Helm et al. (2014). The difference with respect to the processing by ESA mainly consisted of a refined determination of the interferometric phase and of the application of a threshold retracker.

2.3 Accuracy and precision

The accuracy of RA-derived ice surface elevation measurements has been assessed previously by a crossover comparison with independent validation data such as the ICESat laser observations (Brenner et al.2007), airborne lidar (Nilsson et al.2016) and ground-based GNSS profiles (Schröder et al.2017). Besides the offset due to snowpack penetration and instrumental calibration over flat terrain, these assessments revealed that, with increasingly rough surface topography, the RA measurements show systematically higher elevations than the validation data. These topography-related offsets can be explained by the fact that for surfaces that undulate within the ∼20km beam-limited footprint, the radar measurements tend to refer to local topographic maxima (the POCA), while the validation data from ground-based GNSS profiles or ICESat-based profiles represent the full topography. The standard deviation of differences between RA data and validation data contains information about the measurement noise but is additionally influenced by the significantly different sampling of a rough surface as well. While over flat terrain, this standard deviation is below 50 cm for most satellite altimeter data sets, it can reach 10 m and more in coastal regions. However, both types of error relate to the different sampling of topography of the respective observation techniques. An elevation change, detected from within the same technique, is not influenced by these effects. Hence, with respect to elevation changes, not the accuracy but the precision (i.e., the repeatability) has to be considered.

This precision can be studied using intra-mission crossovers between ascending and descending profiles. Here, the precision of a single measurement is obtained by σH=|ΔH|/2 as two profiles contribute to this difference. To reduce the influence of significant real surface elevation changes between the two passes, we only consider crossovers with a time difference of less than 31 days. In stronger inclined topography, the precision of the slope correction dominates the measurement error (Bamber1994). Hence, to provide meaningful results, the surface slope needs to be taken into consideration. We calculate the slope from the CryoSat-2 DEM (Helm et al.2014). The absence of slope-related effects on flat terrain allows for the study of the influence of the retracker (denoted as noise here). With increasing slope, the additional error due to topographic effects can be identified.

A comparison of the crossover errors of our reprocessed data and of the standard products shows significant improvements achieved by our reprocessing. Figure 2 shows this comparison for Envisat (similar plots for each data set can be found in Fig. S1), binned into groups of 0.05 of specific surface slope. The results for a flat topography show that a 10 % threshold provides the highest precision, which confirms the findings of Davis (1997). For higher slopes, we see that our refined slope correction also contributed to a major improvement. A constant noise level σnoise and a quadratic, slope-related term σslope have been fitted to the data according to σH=σnoise+σslopes2, where s is in the unit of degrees. The results in Table 1 show that for each of the PLRA data sets of ERS-1, ERS-2 and Envisat, the measurement noise could be reduced by more than 50 % compared to the ESA product which uses the functional fit retracker ICE-2 (Legrésy and Rémy1997). With respect to the CryoSat-2 standard retracker (Wingham et al.2006a), the improvement is even larger. Improvements are also significant for the slope-related component. For the example of Envisat and a slope of 1, the slope-related component is 1.03 m for the ESA product and only 0.37 m for the reprocessed data. The advanced interferometric processing of the SARIn data achieved similar improvements. For the two early missions Seasat and Geosat, the crossover error of our reprocessed profiles is similar to that of the original data set from GSFC. However, the number of crossover points is significantly increased, especially for Geosat (see Fig. S1). This means that our reprocessing obtained reliable data where the GSFC processor rejected the measurements.

Figure 2Precision of different processing versions of Envisat measurements from near-time (<31 days) crossovers, binned against slope. Red curve: ESA version with ICE-2 retracker and relocated by mean surface slope. Light, medium and dark blue curves: data reprocessed in this study with 50 %-, 20 %- and 10 %-threshold retracker, relocated using the refined method. Vertical bars: number of crossovers for the ESA (red) and our 10 % threshold retracked data (blue).


Table 1Noise level and slope-related component (s in degrees) of the measurement precision, fitted to near-time crossovers (unit: m) of the data from the respective data center and our reprocessed data (with a 10 % threshold retracker applied).

Note that the slope-dependent component is weakly determined for data sets with poor tracking in rugged terrain such as Seasat, Geosat or the ERS ocean mode and for the LRM mode of CryoSat-2.

Download Print Version | Download XLSX

In addition to measurement noise, reflected in the crossover differences, a consistent pattern of offsets between ascending and descending tracks has been observed previously (A–D bias; Legrésy et al.1999; Arthern et al.2001). Legrésy et al. (1999) interpret this pattern as an effect of the interaction of the linearly polarized radar signal with wind-induced surface structures, while Arthern et al. (2001) attribute the differences to anisotropy within the snowpack. Helm et al. (2014) showed that a low threshold retracker significantly reduces the A–D bias. We observe a similar major reduction (from ±1m in some regions for a functional fit retracker to ±15cm when using a 10 % threshold; see Fig. S2). The remaining bias is not larger, in its order of magnitude, than the respective noise. Moreover, near the ice sheet margins, the determination of meaningful A–D biases is complicated by the broad statistical distribution of A–D differences and the difficulty of discriminating outliers. We therefore do not apply a systematic A–D bias as a correction but rather include its effect in the uncertainty estimate of our final result.

3 Multi-mission SEC time series
Back to toptop

3.1 Repeat-track parameter fit

We obtain elevation time series following the repeat-track approach, similar to Legrésy et al. (2006) and Flament and Rémy (2012). As the orbits of the missions used here have different repeat-track patterns, instead of along-track boxes we perform our fit on a regular grid with 1 km spacing (as in Helm et al.2014). For each grid cell we analyze all elevation measurements hi within a radius of 1 km around the grid cell center. This size seems reasonable as for a usual along-track spacing of about 350 m for PLRA (Rémy and Parouty2009), each track will have up to five measurements within the radius. Due to the size of the pulse-limited footprint, a smaller search radius would contain only PLRA measurements with very redundant topographic information and thus would not be suitable to fit a reliable correction for the topography. As specified in Eq. (1), the parameters contain a linear trend (dh∕dt), a planar topography (a0, a1, a2) and a regression coefficient (dBS) for the anomaly of backscattered power (bsi-bs) to account for variations in the penetration depth of the radar signal.

For a single mission, the parameters are adjusted according to the following model:


Here, ti denotes the time of the observation. The reference epoch t0 is set to 09/2010. xi and yi are the polar stereographic coordinates of the measurement location, reduced by the coordinates of the cell's center. The residual resi describes the misfit between the observation and the estimated parameters.

To account for varying penetration depths due to variations in the electromagnetic properties of the ice sheet surface, different approaches exist. Wingham et al. (1998), Davis and Ferguson (2004), McMillan et al. (2014) and Zwally et al. (2015) apply a linear regression using the backscattered power. Flament and Rémy (2012), Michel et al. (2014) and Simonsen and Sørensen (2017) use two additional waveform shape parameters, obtained from functional fit retrackers. Nilsson et al. (2016) showed that a low threshold retracker mitigates the need for a complex waveform shape correction. Hence, we decided to use a solely backscatter-related correction.

Besides the parameters in Eq. (1), McMillan et al. (2014) and Simonsen and Sørensen (2017) estimate an additional orbit-direction-related parameter to account for A–D biases. In Sect. 2.3 we showed that these biases are significantly reduced due to reprocessing with a low threshold retracker. A further reduction of possible remaining artifacts of A–D biases is achieved by smoothing in the merging step in Sect. 3.3.3. The weighted averaging of the results over a diameter of 60 km leads to a balanced ratio of ascending and descending tracks. Our choices concerning the correction for local topography, time-variable penetration effects and A–D biases were guided by the principle of preferring the simplest viable model in order to keep the number of parameters small compared to the number of observations.

In contrast to this single mission approach, here we perform a combined processing of all data from different missions and even different altimeter techniques. Thus, some of the parameters may vary between the data sets. To allow for offsets between the missions, the elevation at the cell center a0 is fitted for each mission individually. The same applies to dBS, which might relate to specific characteristics of a mission as well. For Seasat, covering less than 100 days, this parameter is not estimated, as we assume that during the mission lifetime no significant changes occurred. For ICESat, dBS is not estimated either, as signal penetration is negligible for the laser measurements.

Between different observation techniques (i.e., PLRA, SARIn and laser altimetry), the effective surface slope may also differ. Considering the specific footprint sizes and shapes, the topography is sampled in a completely different way as illustrated in Fig. 3. While PLRA refers to the closest location anywhere within the ∼20km beam-limited footprint (i.e., the POCA), CryoSat's SARIn measurement can be attributed within the narrow Doppler stripe in the cross-track direction. For ICESat the ∼70m laser spot allows a much better sampling of local depressions. Hence, the slope parameters a1 and a2 are estimated for each of the techniques independently.

Considering these sensor-specific differences, the model for the least squares adjustment in Eq. (1) is extended for multi-mission processing to


where M(i) and T(i) denote to which mission or technique the measurement hi belongs.

Figure 3Illustration of the technique-dependent topographic sampling. The laser (red) measures the surface elevation in the nadir of the instrument, while for radar altimetry (blue) the first return signal originates from the POCA (marked by the blue point). Hence, planar surface approximations to the measured heights (dashed lines) as in Eq. (2) are intrinsically different for the different techniques.


We define a priori weights for the measurements hi based on the precision of the respective mission and mode from crossover analysis (Table 1) and depending on the surface slope at the measurement location. This means that in regions with a more distinctive topography, ICESat measurements (with a comparatively low slope-dependent error component) will obtain stronger weights, compared to PLRA as Envisat. Over regions of flat topography, such as the interior of East Antarctica, the weights between PLRA and ICESat are comparable.

In order to remove outliers from the data and the results we apply different outlier filters. After the multi-mission fit, we screen the standardized residuals (Baarda1968) to exclude any resi that exceed 5 times its a posteriori uncertainty. We iteratively repeat the parameter fit until no more outliers are found. Furthermore, in order to exclude remaining unrealistic results from further processing, we filter our repeat-track cells and reject any results where we obtain an absolute elevation change rate |dh/dt| which is larger than 20 m yr−1 or where the standard deviation of this rate is higher than 0.5 m yr−1.

3.2 Single-mission time series

After fitting all parameters according to the multi-mission model (Eq. 2), we regain elevation time series by recombining the parameters a0 and dh∕dt with monthly averages of the residuals (res). For each month j and each mission M, the time series are constructed as


This recombination of parameters from Eq. (2) and averages of residuals does not include the parameters of topography slope and backscatter regression. Hence, each time series of hj,M relates to the cell center and is corrected for time-variable penetration effects. Due to the reference elevation a0,M, which may also contain the inter-mission offset, the penetration depth and a component of the topography sampling within the cell, this results in individual time series for each single mission. A schematic illustration of the results of this step is given in Fig. 5a. The temporal resolution of these time series is defined by using monthly averages of the residuals. With typical repeat-cycle periods of 35 days or more, these res typically represent the anomalies of a single satellite pass towards all parameters including the linear rate of elevation change. The standard deviation of the residuals in these monthly averages is used as uncertainty measure for hj,M (see Sect. S3.1 for further details).

3.3 Combination of the single-mission time series

In order to merge data from different missions into a joint time series, inter-mission offsets have to be determined and eliminated. In the ERS reprocessing project (Brockley et al.2017), mean offsets between the ERS missions and Envisat have been determined and applied to the elevation data. However, for ice sheet studies inter-satellite offsets are found to be regionally varying (Zwally et al.2005; Thomas et al.2008; Khvorostovsky2012). When merging data from different observation techniques (PLRA, SARIn and laser) the calibration gets even more challenging. We chose an approach in different steps which is depicted in Figs. 4 and 5. The following section gives an overview and explains the different steps to merge the single mission time series. A detailed description of the parameters used in each step can be found in the Supplement.

3.3.1 Merging PLRA time series

In the first step, we merge the PLRA time series. For these missions the topographic sampling by the instruments is similar and thus the offsets are valid over larger regions. For overlapping missions (ERS-1, ERS-2, Envisat, CryoSat-2 LRM) the offsets are calculated from simultaneous epochs (blue area in Fig. 5b), as performed by Wingham et al. (1998) or Paolo et al. (2016). Smoothed grids of these offsets are generated, summed up if necessary to make all data sets comparable with Envisat (see Fig. S4) and applied to the respective missions. For the ERS missions, we find significant differences in the offsets for ice and ocean mode; hence, we determine separate offsets for each mode. Comparing our maps with similar maps of offsets between ERS-2 (ice mode) and Envisat shown by Frappart et al. (2016) reveals that the spatial pattern agrees very well, but we find significantly smaller amplitudes. We interpret this as a reduced influence of volume scattering due to our low retracking threshold. In accordance with Zwally et al. (2005), we did not find an appropriate functional relationship between the offset and the waveform parameters.

Figure 4Schematic diagram of the processing steps from the combined repeat-track parameter fit over single-mission time series towards a combined multi-mission time series.


Figure 5Schematic illustration of the combination of the missions. (a) Single-mission time series of PLRA missions (blue and cyan), CryoSat-2 in SARIn mode (green) and the laser altimetry measurements of ICESat (red) with inter-mission offsets. (b) Offsets between the PLRA data are determined from overlapping epochs (blue area) or trend-corrected elevation differences (according to Eq. 2) where dh∕dt is sufficiently stable. (c) The specific offset between PLRA, SARIn and laser data depends on the sampling of the topography within each single cell. These different techniques are aligned by reducing each elevation time series by the specific elevation at the reference epoch tref. Due to possible nonlinear surface elevation changes, this reference elevation is obtained from an 8-year interval only (gray area). (d) The combined multi-mission time series contains SECs with respect to tref.


To calibrate Geosat and Seasat, a gap of several years without observations has to be bridged. As depicted by the dashed blue lines in Fig. 5b, we do this using the trend-corrected reference elevations a0,M from the joint fit in Eq. (2). This, however, can only be done if the rate is sufficiently stable over the whole period. Therefore, we use two criteria. First, we check the standard deviation of the fit of dh∕dt. This σdh∕dt indicates the consistency of the observations towards a linear rate during the observational period. However, anomalies during the temporal gaps between the missions (i.e., 1978–1985 and 1989–1992) cannot be detected in this way. Therefore, we further utilize a firn densification model (FDM; Ligtenberg et al.2011; van Wessem et al.2018). This model describes the anomalies in elevation due to atmospheric processes against the long-term mean. The root mean square (rms) of the FDM time series is hence a good measure for the magnitude of the nonlinear variations in surface elevation. Consequently, only cells where σdh/dt<1 cm yr−1 and rmsFDM<20cm, indicating a highly linear rate, are used to calibrate the two historic missions. Maps of the offsets with respect to Envisat are shown in Fig. S5. The FDM criterion is not able to detect changes in ice dynamics. However, as regions where both stability criteria are fulfilled are mainly found on the plateau where flow velocity is below 30 m yr−1 (Rignot et al.2017), we expect no significant nonlinear elevation changes due to ice flow. The mean of the offsets over all cells amounts to −0.86m for Seasat and −0.73m for Geosat. The corresponding standard deviations of 0.85 and 0.61 m are mainly a result of the regional pattern of the offsets. The true offsets are likely to have spatial variations. However, we are not able to distinguish spatial variations of the offset from residual effects of temporal height variations in the regions meeting the stability criterion. In the regions not meeting this criterion, we are not able to estimate the spatial variations of the offset at all. Therefore, our final estimate of the offset, applied to the measurements, is a constant, calculated as the average offset over the regions meeting the stability criterion. The spatial variability not accounted for by the applied offset is included, instead, in the assessed uncertainty. Our bias between Seasat and Envisat (-0.86±0.85m) agrees within uncertainties with the ocean-based bias of −0.77m used by Fricker and Padman (2012). However, we prefer the offset determined over the ice sheet because this kind of offset may depend on the reflecting medium (see Sect. S3.2.2 for a more detailed discussion).

With the help of these offsets, all PLRA missions were corrected towards the chosen reference mission Envisat. Uncertainty estimates of the offsets are applied to the respective time series to account for the additional uncertainty. Hence, the PLRA time series are combined (blue in Fig. 5c with additional CryoSat-2 LRM mode where available). At epochs when more than one data set exists, we apply weighted averaging using the uncertainty estimates.

3.3.2 Technique-specific surface elevation changes

In contrast to the PLRA data in the previous step, when merging data from different observation techniques such as CryoSat's SARIn mode, ICESat's laser observations and PLRA, the different sampling of topography also has to be considered. As noted in Sect. 3.1, this might lead to completely different surfaces fitted to each type of elevation measurement and the time series need to be calibrated for each cell individually. However, not all cells have valid observations of each data set. Therefore, instead of calibrating the techniques against each other, we reduce each time series by their elevation at a common reference epoch and hence obtain time series of surface elevation changes (SEC) with respect to this reference epoch instead of absolute elevation time series. This step eliminates offsets due to differences in firn penetration or due to the system calibration between the techniques as well.

We chose September 2010 (09/2010) as the reference epoch. This epoch is covered by the observational periods of PLRA and CryoSat SARIn and also is exactly 1 year after the last observations of ICESat, which reduces the influence of an annual cycle. As discussed in Sect. 3.3.1, nonlinear elevation changes will adulterate a0 from Eq. (2), obtained over the full time span. Therefore, we applied another linear fit to a limited time interval of 8 years only (September 2006–September 2014, gray area in Fig. 5c). We subtract the variation of the FDM over this period to account for short-term variations. The limited time interval reduces the influence of changes in ice dynamics. We estimate the individual reference elevations a0,T for each technique T and a joint dh∕dt. After subtracting the technique-specific reference elevations a0,T from the respective time series, they all refer to 09/2010 and can be combined.

Figure 6Five example snapshots of the resulting combined surface elevation time series (a) and their corresponding uncertainty (b). The height differences refer to our reference epoch 09/2010.


3.3.3 Merging different techniques

We perform the final combination of the techniques using a weighted spatiotemporal averaging with 10 km σ Gaussian weights in the spatial domain (up to a radius of 3σ=30km) and over 3 epochs (i.e., including the two consecutive epochs) in the temporal domain. Hence, we obtain grids of surface elevation change (SEC) with respect to 09/2010 for each month observed. Due to the smoothing of the weighting function, we reduce our spatial SEC grid resolution to 10 km× 10 km. The respective uncertainties are calculated according to the error propagation. To avoid extrapolation and to limit this merging step to the observed area only, we calculate a value for an epoch in the 10 km× 10 km grid cells only if we have data within 20 km around the cell center (which is about the size of a beam-limited radar footprint). The five examples in Fig. 6 demonstrate the spatiotemporal coverage of the resulting SEC grids at different epochs. The corresponding uncertainty estimates, given in Fig. 6b (further details in the Supplement), reach values of 1 m and more. Besides the measurement noise and the uncertainties of the offsets, these uncertainty estimates contain a further component from the weighted averaging. In regions with large variations in Δh over relatively short spatial scales (such as at fast-flowing outlet glaciers), such variations can add a significant contribution to σΔh. As the magnitude of Δh grows with the temporal distance to the reference epoch, the largest contributions to σΔh can be expected for the earliest epochs. This also explains why the epoch 09/2008 provides the lowest uncertainty estimate in these examples, even lower than the CryoSat-2-based epoch 09/2017.

Figure 7Validation with elevation differences observed by kinematic GNSS between 2001 and 2015 (a, b, c) and Operation IceBridge between 2002 and 2016 (d, e, f). Differences between elevation changes observed by the validation data and altimetry are shown in the maps (a, d, color scale in panels b and e). Median and MAD of these differences, binned by different surface slope, are shown in the center (b, e). The right diagrams (c, f) show the comparison of these differences with the respective uncertainty estimate, obtained from both data sets. The point density is plotted from yellow to blue and the black dots show the root mean square, binned against the estimated uncertainty.


4 Comparison of SEC with independent data
Back to toptop

4.1 In situ and airborne observations

To validate our results, we used inter-profile crossover differences of 19 kinematic GNSS profiles (available for download from Schröder et al.2016) and elevation differences from Operation IceBridge (OIB ATM L4; Studinger2014). The ground-based GNSS profiles were observed between 2001 and 2015 on traverse vehicles of the Russian Antarctic Expedition and most of them covered more than 1000 km. The accuracy of these profiles has been determined in Schröder et al. (2017) to 4–9 cm. One profile (K08C) has not been used due to poorly determined antenna height offsets. For each crossover difference between kinematic profiles from different years, we compare the differences of the corresponding altimetric SEC epochs in this location (δΔh=ΔhKIN-ΔhALT). The same analysis has been performed with the elevation changes obtained from differences of measurements of the scanning laser altimeter (Airborne Topographic Mapper, ATM) of OIB. As described by Studinger (2014), the Level 4 Δh product is obtained by comparing planes fitted to the laser scanner point clouds. The flights, carried out between 2002 and 2016, were strongly concentrated along the outlet glaciers of West Antarctica and the Antarctic Peninsula. Hence, they cover much more rugged terrain, which is more challenging for satellite altimetry. Over the tributaries of the Amundsen Sea glaciers and along the polar gap of ICESat, some repeated measurements have also been performed over flat terrain. The accuracy of these airborne measurements has been validated, e.g., near summit station in Greenland. Brunt et al. (2017) used ground-based GNSS profiles of snowmobiles for this task and obtained offsets on the order of only a few centimeters and standard deviations between 4 and 9 cm.

Figure 7a and d show the results of our validation (more detailed maps for several regions at Fig. S6). A satellite calibration error would lead to systematic biases between the observed elevation differences if ΔhALT is obtained from data of two different missions. However, such biases may also be caused by systematic errors in the validation data. Furthermore, in contrast to the calibration data, the RA measurements may systematically miss some of the most rapid changes if those are located in local depressions (Thomas et al.2008). With an overall median difference of 6±10cm for the GNSS profiles and -9±42cm for OIB, however, the observed elevation changes show only moderate systematic effects and agree within their error bars. The median absolute deviation (MAD) for different specific surface slopes (Fig. 7b and e) reveals the influence of topography in this validation. The GNSS profiles show only a very small increase of this variation with slope. The IceBridge data cover the margins of many West Antarctic glaciers, where elevation changes differ over relatively short distances. Hence, it is not surprising that we see a significantly larger spread of δΔh at higher slopes here. However, in the less inclined regions, the MAD of the differences is still at a level of 25 cm, which is significantly larger than in the comparison with the GNSS profiles. This large spread for regions with low slopes originates mainly from the tributaries of Pine Island Glacier, where many campaigns of OIB are focused (see Figs. 7d and S6d for details). While still relatively flat, the surface elevation in this area is comparatively low, which leads to a stronger influence of precipitation (Fyke et al.2017). This induces higher short-term variations in surface elevation, which might explain the higher differences between the IceBridge results and our 3-month temporally smoothed altimetry data. In contrast, the differences around the South Pole or in Queen Elizabeth Land (see also Fig. S6c) are significantly smaller. For the 2016 campaign of OIB, Brunt et al. (2018) furthermore report a spurious elevation variation of 10 to 15 cm across the wide-scan ATM swath, which indicates a bias in the instrumental tilt angle. This could explain the systematic differences along the 88 S circle profiles where this campaign is involved.

The observed δΔh can further be used to evaluate the uncertainty estimates. In Fig. 7c and f, the uncertainty estimates of the four contributing data sets are combined and compared to the observed differences. The comparison with both validation data sets supports the conclusion that the uncertainty estimates are reasonable. For ΔhALT we expect higher errors in coastal regions due to the increased uncertainty of the topographic correction in radar altimetry. A similar relation to topography is expected for ΔhOIB due to the plane fit to the ATM point cloud, but surface roughness and crevassing also play an important role here. In contrast, the errors of the GNSS-derived ΔhKIN are almost independent of topography. Instead, ΔhKIN tends to be more uncertain on the plateau, where the soft snow causes large variations of the subsidence of the vehicles into the upper firn layers. The relatively low differences in δΔh, even in regions that imply a higher uncertainty, are likely just incidental for the small sample of validation data along the GNSS profiles.

In conclusion, this validation shows that remaining systematic biases (originating from satellite altimetry or the validation data) are less than a decimeter in the observed regions and that our uncertainty estimate is realistic. However, only altimetric SEC within the interval 2001–2016 can be validated in this way. For the earlier missions, no spatially extensive high-precision in situ data are publicly available.

Figure 8(a) Correlation coefficient between the SEC anomalies of the altimetry grids and the FDM over 1992–2016 after detrending. (b) Average magnitude of anomalies of the altimetry time series. (c) Correlation coefficient plotted against the nonlinear SEC. The point density is color coded from yellow to blue. The black dots show the binned mean values.


4.2 Firn model

Another data set, which covers almost the identical spatial and temporal range as the altimetric data, is the IMAU Firn Densification Model (FDM; Ligtenberg et al.2011), forced at the upper boundary by the accumulation and temperature of the Regional Atmospheric Climate Model, version 2.3p2 (van Wessem et al.2018). The IMAU-FDM has been updated to the period 1979–2016, modeling the firn properties and the related surface elevation changes on a 27 km grid. However, as the FDM contains elevation anomalies only, any long-term elevation trend over 1979–2016, e.g., due to changes in precipitation on longer timescales (for example those observed in some regions of West Antarctica; Thomas et al.2015) would not be included in the model. Furthermore, due to the nature of the model, it cannot give information about ice dynamic thinning and thickening. Hence, to compare the FDM and the SEC from altimetry, we first remove a linear trend from both data sets. This is performed for the period 1992–2016 (depicted in Fig. S7). The trends are only calculated from epochs where both data sets have data; i.e., in the polar gap this comparison is limited to 2003–2016 or 2010–2016, depending on the first altimetry mission providing data here. After the detrending, the anomalies are used to calculate correlation coefficients for each cell, depicted in Fig. 8a. Figure 8b shows the average magnitude of the seasonal and interannual variations (nonlinear SEC), calculated as the rms of the anomalies from the altimetry data. Comparing the two maps shows that the correlation is around 0.5 or higher, except in regions where the magnitude of the anomalies is small (i.e., where the signal-to-noise ratio of the altimetric data is low) and where large accelerations in ice velocity are observed (such as near the grounding zone of Pine Island Glacier). The relationship between the correlation coefficient and the magnitude of the nonlinear SEC is depicted in Fig. 8c, where we see that for the vast majority of cells the correlation is positive. For anomalies with a nonlinear SEC >0.5m, the average correlation is between 0.3 and 0.6.

Figure 9Multi-mission surface elevation change from the combined SEC time series over different time intervals. (a, b) The long-term surface elevation change between 1978 and 2017 and 1992 and 2017 for the covered area. (c–j) Elevation change over consecutive time intervals reveals the interannual variability. Thin lines mark the drainage basin outlines, denoted in panel (a). Bold letters in boxes in panel (b) denote areas mentioned in the text and in Fig. 10.


Anomalies against the simultaneously observed long-term trend (1992–2016) can also be computed for earlier epochs. Assuming no significant changes in ice dynamics here, these anomalies allow for a comparison of Geosat and Seasat with the FDM. The median difference between the anomalies according to Geosat and the anomalies according to the FDM amounts to 0.12±0.21m (see Fig. S8). Considering that this difference is very sensitive to extrapolating the long-term trends, this is a remarkable agreement. With a median of 0.26±0.32m, the difference between anomalies from Seasat and from the FDM is larger, but this comparison is also more vulnerable to potential errors due to the extrapolation. As the FDM starts in 1979 while Seasat operated in 1978, we compare the Seasat data with the FDM anomalies from the respective months of 1979, which might impose additional differences. Finally, the FDM model has its own inherent errors and uncertainties. Therefore, only part of the differences originate from errors in the altimetry results.

5 Results
Back to toptop

5.1 Surface elevation changes

The average rates of elevation change over different time intervals of our multi-mission time series are shown in Fig. 9. To calculate these rates, we first averaged the data over the first year and the last year of the interval to reduce the noise, then subtracted the respective averages from each other and finally divided these differences by their time difference in years. If one of the years does not cover the full annual cycle, we calculate the average only from the months covered in both years (July–October for 1978–2017, April–December for 1992–2017). We calculate the SEC rate from epoch differences instead of fitting a rate to all epochs because the first observations at specific latitudes start in different years, the observations have different precisions and the large gap between 1978 and 1985 is not covered by observations at all. These three points would lead to a bias towards the later epochs in a fit, so that the rates would not be representative for the true average elevation change over the full interval.

The long-term elevation changes over 25 years (Fig. 9b) show the well known thinning in the Amundsen Sea Embayment and at Totten Glacier, as well as the thickening of Kamb Ice Stream (see, e.g., Wingham et al.2006b; Flament and Rémy2012; Helm et al.2014). In contrast, 60 % of East Antarctica north of 81.5 S shows surface elevation changes of less than ±1 cm yr−1. Several coastal regions of the EAIS, however, show significant elevation changes. Totten Glacier (T in Fig. 9b) is thinning at an average rate of 72±18 cm yr−1 at the grounding line (cf. Fig. 10b). Several smaller glaciers in Wilkes Land also show a persistent thinning. We observe SEC rates of -26±10 cm yr−1 at Denman Glacier (D), -41±19 cm yr−1 at Frost Glacier (F) and -33±12 cm yr−1 near Cook Ice Shelf (C). Rignot (2006) showed that the flow velocity of these glaciers, which are grounded well below sea level, was above the balance velocity for many years. Miles et al. (2018) analyzed satellite images since 1973 and found that the flow velocity of Cook Ice Shelf glaciers has significantly accelerated since then. In contrast, the western sector of the EAIS (Coats Land, Dronning Maud Land and Enderby Land; basins J”–B) shows thickening over the last 25 years at rates of up to a decimeter per year.

Figure 10Multi-mission SEC time series in four selected regions (a) Pine Island Glacier, (b) Totten Glacier, (c) Shirase Glacier in Dronning Maud Land and (d) Lake Vostok (marked by P, T, S and V in Fig. 9b). The time series of points B, C and D are shifted along Δh for better visibility and the one σ uncertainty range is displayed in black. The maps on the left show the elevation change rate between 1992 and 2017 as in Fig. 9b (but in a different color scale).


Comparing the long-term elevation changes over 40 years (Fig. 9a) with those over 25 years shows the limitations of the early observations, but also the additional information they provide. There were relatively few successful observations at the margins. However, for Totten and Denman glaciers, the 40-year rates at a distance of approximately 100 km inland from the grounding line are similar to the rates over the 1992–2017 interval, which indicates a persistent rate of thinning. Another benefit of our merged time series is that they allow for the calculation of rates over any subinterval, independent of mission periods as demonstrated in Fig. 9c–j. For most of the coastal regions of the AIS, these rates over different intervals reveal that there is significant interannual variation. Such large-scale fluctuations in elevation change have been previously reported by Horwath et al. (2012) and Mémin et al. (2015) for the Envisat period. Our combined multi-mission time series now allow a detailed analysis of such signals on a temporal scale of up to 40 years.

Four examples for elevation change time series in the resulting multi-mission SEC grids are shown in Fig. 10 (coordinates in Table S2). Pine Island Glacier (PIG) is located in the Amundsen Sea Embayment, which is responsible for the largest mass losses of the Antarctic Ice Sheet (e.g., McMillan et al.2014). In East Antarctica, the largest thinning rates are observed at Totten Glacier. The region of Dronning Maud Land and Enderby Land in East Antarctica has been chosen as an example for interannual variation. Here, Boening et al. (2012) reported two extreme accumulation events in 2009 and 2011, which led to significant mass anomalies. We chose a profile at Shirase Glacier as an example for this region. In contrast to the previous locations, a very stable surface elevation has been reported for Lake Vostok (e.g., Richter et al.2014). This stability, however, has been a controversial case recently (Zwally et al.2015; Scambos and Shuman2016; Richter et al.2016). Therefore, our results in this region shall add further evidence that pinpoints the changes there.

For Pine Island Glacier (Fig. 10a), we observe a continuous thinning over the whole observational period since 1992 (Seasat and Geosat measurements did not cover this region). Close to the grounding line (point D) the surface elevation has decreased by -45.8±7.8m since 1992, which means an average SEC rate of -1.80±0.31 m yr−1. The time series reveals that this thinning was not constant over time, but accelerated significantly around 2006. The mean rate at D over 1992–2006 of -1.32±0.66 m yr−1 increased to -4.17±1.67 m yr−1 over 2007–2010. After 2010, the thinning rates near the grounding line decelerate again and for the period 2013–2017, the rate at D of -1.31±0.80 m yr−1 is very close to the rate preceding the acceleration. Also, at greater distances from the grounding line (B at 80 km, A at 130 km) we observe an acceleration of the prevailing rates around 2006 (-0.44±0.15 m yr−1 over 1992–2006, -1.20±0.10 m yr−1 over 2006–2017 at A). In contrast to the points near the grounding line, further inland the thinning has not decelerated so far and is still at a high level. Hence, for the most recent period (2013–2017) the elevation at all points along the 130 km of the main flow line is decreasing at very similar rates. A similar acceleration of the elevation change rate near the grounding line, followed by slowdown, is observed by Konrad et al. (2016). The onset of this acceleration coincides with the detaching of the ice shelf from a pinning point (Rignot et al.2014). For the time after 2009, Joughin et al. (2016) report relatively little grounding line migration, resulting in a leveling off of the ice flow velocity. This agrees with our observed slowdown of elevation changes.

For Totten Glacier in East Antarctica (Fig. 10b), we observe a clear negative SEC. This has been previously reported by several authors (e.g., Pritchard et al.2009; Flament and Rémy2012; Zwally et al.2015), but our data provide an unprecedented time span and temporal resolution, allowing for the analysis of the evolution of the elevation changes on a monthly scale over up to 40 years. At the grounding line (point D), Totten Glacier thinned by 31.8±7.7m between 1987 and 2017, which results in an average SEC rate of -1.03±0.25 m yr−1. Seasat could not provide successful observations at the grounding line, but the time series for point C (around 60 km inland) with a rate of -0.38±0.10 m yr−1 between 1978 and 2017 and for point B (150 km) with a rate of -0.11±0.04 m yr−1 indicate that this thinning preceded the epoch of Geosat. At point A in a distance of 280 km, we find no significant elevation change (0.01±0.03 m yr−1 for 1978–2017). The temporal resolution of thes