- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications

Journal cover
Journal topic
**The Cryosphere**
An interactive open-access journal of the European Geosciences Union

Journal topic

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications

**Research article**
05 Feb 2019

**Research article** | 05 Feb 2019

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

^{1}Technische Universität Dresden, Institut für Planetare Geodäsie, Dresden, Germany^{2}Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany^{3}Institute for Marine and Atmospheric Research Utrecht, Utrecht University, Utrecht, the Netherlands

^{1}Technische Universität Dresden, Institut für Planetare Geodäsie, Dresden, Germany^{2}Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany^{3}Institute for Marine and Atmospheric Research Utrecht, Utrecht University, Utrecht, the Netherlands

**Correspondence**: Ludwig Schröder (ludwig.schroeder@tu-dresden.de)

**Correspondence**: Ludwig Schröder (ludwig.schroeder@tu-dresden.de)

Abstract

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 $-\mathrm{85}\pm \mathrm{16}$ Gt yr^{−1} between 1992 and
2017, which accelerated to $-\mathrm{137}\pm \mathrm{25}$ Gt yr^{−1} after 2010.

How to cite

Back to top
top
How to cite.

Schröder, L., Horwath, M., Dietrich, R., Helm, V., van den Broeke, M. R., and Ligtenberg, S. R. M.: Four decades of Antarctic surface elevation changes from multi-mission satellite altimetry, The Cryosphere, 13, 427–449, https://doi.org/10.5194/tc-13-427-2019, 2019.

1 Introduction

Back to toptop
Satellite altimetry is fundamental for detecting and understanding changes in
the Antarctic Ice Sheet (AIS; Rémy and Parouty, 2009; 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émy, 2012). 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.

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 https://doi.pangaea.de/10.1594/PANGAEA.897390
(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*σ*=20 km. 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 Data

Back to toptop
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).

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 (Bamber, 1994). 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 (Bamber, 1994), 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
(Bamber, 1994). 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 *σ*=1 km 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 (Bamber, 1994). Functional fit approaches (e.g., Martin et al., 1983; Davis, 1992; 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 (Davis, 1997; 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 (Bamber, 1994). 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.

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 ∼20 km 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 ${\mathit{\sigma}}_{H}=\left|\mathrm{\Delta}H\right|/\sqrt{\mathrm{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 (Bamber, 1994). 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
${\mathit{\sigma}}_{H}={\mathit{\sigma}}_{\mathrm{noise}}+{\mathit{\sigma}}_{\mathrm{slope}}\cdot {s}^{\mathrm{2}}$, 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émy, 1997).
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.

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.

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 ±1 m in some regions for a functional fit retracker to ±15 cm 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
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 *h*_{i} 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 Parouty, 2009), 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 (d*h*∕d*t*), a planar
topography (*a*_{0}, *a*_{1}, *a*_{2}) and a regression coefficient (dBS) for the
anomaly of backscattered power (bs${}_{i}-\stackrel{\mathrm{\u203e}}{\mathrm{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:

$$\begin{array}{ll}{\displaystyle}{h}_{i}=& {\displaystyle}\mathrm{d}h/\mathrm{d}t({t}_{i}-{t}_{\mathrm{0}})+\\ {\displaystyle}& {\displaystyle}{a}_{\mathrm{0}}+{a}_{\mathrm{1}}{x}_{i}+{a}_{\mathrm{2}}{y}_{i}+\\ {\displaystyle}& {\displaystyle}\mathrm{dBS}\left({\mathrm{bs}}_{i}-\stackrel{\mathrm{\u203e}}{\mathrm{bs}}\right)+\\ \text{(1)}& {\displaystyle}& {\displaystyle}{\mathrm{res}}_{i}.\end{array}$$

Here, *t*_{i} denotes the time of the observation. The reference epoch *t*_{0} is
set to 09/2010. *x*_{i} and *y*_{i} are the polar stereographic coordinates of
the measurement location, reduced by the coordinates of the cell's center.
The residual res_{i} 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
*a*_{0} 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
∼20 km 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 ∼70 m
laser spot allows a much better sampling of local depressions. Hence, the
slope parameters *a*_{1} and *a*_{2} 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

$$\begin{array}{ll}{\displaystyle}{h}_{i}=& {\displaystyle}\mathrm{d}h/\mathrm{d}t({t}_{i}-{t}_{\mathrm{0}})+\\ {\displaystyle}& {\displaystyle}{a}_{\mathrm{0},M\left(i\right)}+{a}_{\mathrm{1},T\left(i\right)}{x}_{i}+{a}_{\mathrm{2},T\left(i\right)}{y}_{i}+\\ {\displaystyle}& {\displaystyle}{\mathrm{dBS}}_{M\left(i\right)}\left({\mathrm{bs}}_{i}-\stackrel{\mathrm{\u203e}}{{\mathrm{bs}}_{M\left(i\right)}}\right)+\\ \text{(2)}& {\displaystyle}& {\displaystyle}{\mathrm{res}}_{i},\end{array}$$

where *M*(*i*) and *T*(*i*) denote to which mission or technique the measurement
*h*_{i} belongs.

We define a priori weights for the measurements *h*_{i} 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 (Baarda, 1968) to exclude any res_{i} 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
$|\mathrm{d}h/\mathrm{d}t|$ which is larger than 20 m yr^{−1} or where the standard
deviation of this rate is higher than 0.5 m yr^{−1}.

After fitting all parameters according to the multi-mission model
(Eq. 2), we regain elevation time series by recombining the
parameters *a*_{0} and d*h*∕d*t* with monthly averages of the
residuals ($\stackrel{\mathrm{\u203e}}{\mathrm{res}}$). For each month *j* and each mission
*M*, the time series are constructed as

$$\begin{array}{}\text{(3)}& {h}_{j,M}={a}_{\mathrm{0},M}+\mathrm{d}h/\mathrm{d}t({t}_{j}-{t}_{\mathrm{0}})+\stackrel{\mathrm{\u203e}}{{\mathrm{res}}_{j,M}}.\end{array}$$

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 *h*_{j,M} relates to the cell center and is corrected for
time-variable penetration effects. Due to the reference elevation *a*_{0,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
$\stackrel{\mathrm{\u203e}}{\mathrm{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 *h*_{j,M} (see Sect. S3.1 for further
details).

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; Khvorostovsky, 2012). 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.

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.

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 *a*_{0,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 d*h*∕d*t*. 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
${\mathit{\sigma}}_{\mathrm{d}h/\mathrm{d}t}<\mathrm{1}$ cm yr^{−1} and
rms_{FDM}<20 cm, 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.86 m for Seasat and −0.73 m 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 ($-\mathrm{0.86}\pm \mathrm{0.85}$ m)
agrees within uncertainties with the ocean-based bias of −0.77 m
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.

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
*a*_{0} 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 *a*_{0,T} for each
technique *T* and a joint d*h*∕d*t*. After subtracting the
technique-specific reference elevations *a*_{0,T} from the respective time
series, they all refer to 09/2010 and can be combined.

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*σ*=30 km) 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.

4 Comparison of SEC with independent data

Back to toptop
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; Studinger, 2014). 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 ($\mathit{\delta}\mathrm{\Delta}h=\mathrm{\Delta}{h}_{\mathrm{KIN}}-\mathrm{\Delta}{h}_{\mathrm{ALT}}$). 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
Δ*h*_{ALT} 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±10 cm for the GNSS profiles and $-\mathrm{9}\pm \mathrm{42}$ cm 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 Δ*h*_{ALT} 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 Δ*h*_{OIB} 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 Δ*h*_{KIN} are almost
independent of topography. Instead, Δ*h*_{KIN} 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.

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.5 m, the average correlation is between 0.3 and 0.6.

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.21 m (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.32 m, 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
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émy, 2012;
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 $-\mathrm{26}\pm \mathrm{10}$ cm yr^{−1} at Denman Glacier (D), $-\mathrm{41}\pm \mathrm{19}$ cm yr^{−1} at Frost
Glacier (F) and $-\mathrm{33}\pm \mathrm{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.

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 Shuman, 2016; 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 $-\mathrm{45.8}\pm \mathrm{7.8}$ m since 1992,
which means an average SEC rate of $-\mathrm{1.80}\pm \mathrm{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
$-\mathrm{1.32}\pm \mathrm{0.66}$ m yr^{−1} increased to $-\mathrm{4.17}\pm \mathrm{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
$-\mathrm{1.31}\pm \mathrm{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 ($-\mathrm{0.44}\pm \mathrm{0.15}$ m yr^{−1} over 1992–2006,
$-\mathrm{1.20}\pm \mathrm{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émy, 2012; 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.7 m between 1987 and 2017, which results in an average
SEC rate of $-\mathrm{1.03}\pm \mathrm{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
$-\mathrm{0.38}\pm \mathrm{0.10}$ m yr^{−1} between 1978 and 2017 and for point B
(150 km) with a rate of $-\mathrm{0.11}\pm \mathrm{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