Four decades of surface elevation change of the Antarctic Ice Sheet from multi-mission satellite altimetry

We developed an approach for a multi-mission satellite altimetry analysis over the Antarctic Ice Sheet which comprises Seasat, Geosat, ERS-1, ERS-2, Envisat, ICESat and CryoSat-2. In a first step we apply a consistent reprocessing of the radar alitmetry data which improves the measurement precision by up to 50%. We then perform a joint repeat altimetry analysis of all missions. We estimate inter-mission offsets by approaches adapted to the temporal overlap or non-overlap and to the similarity or dissimilarity of involved altimetry techniques. Hence, we obtain monthly grids forming a combined surface 5 elevation change time series. Owing to the early missions Seasat and Geosat, the time series span almost four decades from 07/1978 to 12/2017 over 25% of the ice sheet area (coastal regions of East Antarctica and the Antarctic Peninsula). Since the launch of ERS-1 79% of the ice sheet area is covered by observations. Over this area, we obtain a negative volume trend of -34±5 km3/yr for the more than 25-year period (04/1992–12/2017). These volume losses have significantly accelerated to a rate of -170±11 km3/yr for 2010–2017. Interannual variations significantly impact decadal volume rates which highlights 10 the importance of the long-term time series. Our time series show a high coincidence with modeled cumulated precipitation anomalies and with satellite gravimetry. This supports the interpretation with respect to snowfall anomalies or dynamic thinning. Moreover, the correlation with cumulated precipitation anomalies back to the Seasat and Geosat periods highlights that the inter-mission offsets were successfully corrected and that the early missions add valuable information.


Introduction
Satellite altimetry allows to observe the surface elevation changes of the Antarctic Ice Sheet with unprecedented precision and resolution (Rémy and Parouty, 2009;Shepherd et al., 2012).Different missions concurrently observed the dynamic thinning of several outlet glaciers in West Antarctica and the relative stability of most regions of the East Antarctic Ice Sheet (e.g.Wingham et al., 2006;Flament and Rémy, 2012;Helm et al., 2014).However, when going into detail, some significant differences between the rates obtained from different time intervals become evident.While Wingham et al. (2006) observed negative elevation rates in eastern Dronning Maud Land and Enderby Land (25-60°E) and positive rates in Princess Elizabeth Land (70-100°E) using ERS-1 and ERS-2 (1992ERS-2 ( -2003)), Flament and Rémy (2012) observed a contrary pattern using Envisat (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010).Lenaerts et al. (2013)  Such events can have a strong influence on the calculation of surface elevation change rates.Only an observation time span as long as possible can help to reduce the influence of such events and obtain the long-term elevation trend.
Missions with similar sensor characteristics have successfully been combined by Wingham et al. (2006, ERS-1 and ERS-2) or Li and Davis (2008, ERS-2 and Envisat).However, as demonstrated by Khvorostovsky (2012), the correction of intermission biases over an ice sheet is not trivial.When including missions with different sensor characteristics (as ICESat or CryoSat-2 SARIn), a thorough calibration becomes even more important.
Here, we present an approach to combine the different satellite altimetry missions in order to extend the observation time span as long as possible.We create a combined time series of surface elevation change (SEC) that allows to identify rapid changes associated, e.g., to snowfall events as well as long-term changes as e.g.due to changing ice dynamics over nearly four decades.

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 over 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 ERS-1, ERS-2, Envisat and CryoSat-2 the most recent ESA products were used.For ICESat the final release from the National Snow and Ice Data Center (NSIDC) was employed.The inter-campaign biases between the ICESat laser operation periods were corrected following Schröder et al. (2017).To remove corrupted measurements, we edited each data set in a preprocessing step.
Further information concerning the dataset versions used and details about the flags and thresholds applied in the data editing can be found in the supplementary material.
ERS-1 and ERS-2 measurements were performed in two different modes.With a larger tracking window, the measurements in ice mode are more reliable in coastal areas.However, a remarkable amount of observations has been performed in ocean mode over Antarctica as well (22% for ERS-1, 2% for ERS-2).Therefore, similar to Paolo et al. (2016), we use the data from both modes.However, we treat them as two separate data sets as there is a regionally varying bias between the modes.version with ICE-2 retracker and relocated by mean surface slope, light, medium and dark blue stands for our reprocessed data with 50%-, 20%-and 10%-threshold retracker, respectively, relocated using our refined approach.The bars in the background indicate the number of crossovers for the ESA (red) and our 10% (blue) data.

Accuracy and precision
The accuracy of radar altimetry derived ice surface elevation measurements has been assessed previously by a crossover comparison with independent data such as the ICESat laser observations (Brenner et al., 2007) or ground based GNSS profiles (Schröder et al., 2017).Such assessments revealed biases which are highly related to the topography.Due to the fact that a pulse limited radar altimeter measures always the POCA, which is the local maximum within the ~20 km footprint, but a GNSS-or ICESat-based profile represents the full topography, these differences always include a positive bias over undulating surfaces.
Hence, the radar altimetry profiles tend to get biased into positive direction and the standard deviations reach ten meters and more in distinct topography.Nevertheless, this bias does not influence the detection of elevation changes from a single mission.
Hence, with respect to the detection of elevation changes not the accuracy but the precision has to be considered.
This precision (i.e.repeatability) can be studied using intra-mission crossovers between ascending and descending profiles.
Here, the precision of a single measurement is obtained by dividing the absolute value of the crossover difference between two profiles by √ 2. To rule out significant surface elevation changes between the two passes, we consider only 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 to study the influence of the retracker.With increasing slope, the additional error due to topographic effects can be identified.
A comparison of the crossover errors of our reprocessed data with the respective results of the processing versions from the different data centers shows which significant improvements could be achieved by the reprocessing steps described above (see Fig. 2 for Envisat, similar plots for each dataset can be found in the supplementary material).The results for a flat topography show that a 10% threshold provides the highest precision, confirming the findings of Davis (1997).For each of the PLRA datasets of ERS-1, ERS-2 and Envisat, the crossover error could be reduced by about 50% compared to ESA's ICE-2 fit.With respect to the CryoSat-2 CFI retracker, the improvement is even larger.With increasing influence of the topography, our refined Simonsen and Sørensen (2017) estimate an orbit direction related parameter in their repeat track processing to remove this effect, while Helm et al. (2014) showed, that a low threshold retracker significantly reduces the ascending-descending bias.We observe a similar major reduction (from ±1 m in some regions for a functional fit retracker to less than a decimeter when using a 10% threshold, see Fig. S1 in the supplements).Since, moreover, the final results are smoothed over several kilometers and thus usually contain ascending and descending satellite tracks alike, we conclude that this bias has no discernible effect on our results.
The crossover comparison is not only performed for quality control of our processing chain, we also use the results to adequately set weights when combining measurements from different missions in a certain location in the repeat altimetry fit.
Therefore we bin the single profile crossover errors (|∆h|/ √ 2) into 20 groups of specific surface slope and fit a constant plus a quadratic, slope related term to the bin medians ( weights of the small footprint measurements of ICESat in regions with a more distinctive topography.Instead, over the flat interior of East Antarctica, the weights are quite similar.
3 Multi-mission analysis

Surface elevation change time series
We obtain elevation time series following the repeat altimetry 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, we perform our fit on a regular grid with 1 km spacing (as in Helm et al., 2014), instead of along-track boxes.For each grid cell we analyze all elevation measurements h i within a radius of 1 km around the grid cell center.As specified in Eq. 1, we fit a linear trend (dh/dt), a plane (a 0 , a 1 , a 2 ) and a regression coefficient (dBS) for the anomaly of backscattered power (bs i − bs).The search radius of 1 km 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 5 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.
For a single mission, the parameters are adjusted according to the model 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.
For the combined processing of different missions and different altimeter techniques, some of the parameters may vary between the datasets.Thus, they are estimated individually for the respective mission M (i) or observation technique T (i) of the measurement h i .Hence, the general relation for a combined processing can be expressed as To allow for offsets between the missions, the elevation at the cell center a 0 is fitted for each mission M (i) individually.The same is applied 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 life time no significant changes occurred.The same applies to ICESat, where signal penetration is negligible. 6 The Cryosphere Discuss., https://doi.org/10.Between different observation techniques (i.e.PLRA, SARIn and laser altimetry), also the effective surface slope may 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 very narrow Doppler stripe in across track direction.For ICESat the ~70 m laser spot allows a much better sampling of local depressions too.Hence, the slope parameters a 1 and a 2 are estimated for each of the techniques T (i) independently.This setting allows to fit dh/dt from all missions while still considering the individual peculiarities of each sensor.Also in the weighting of the observations we consider their different characteristics.Using the sensor-specific precision (Tab. 1) we set the weights according to the surface slope at the respective measurements location.
Outliers are removed iteratively in the processing.Therefore, we exclude observations where the standardized residuals exceed a value of 5 and repeat the processing until no more outliers are found.
After fitting all parameters according to Eq. 2, we can regain elevation time series by recombining the parameters with the residuals.We use monthly averages of the residuals, which typically represent the misfit of a single satellite pass towards all respective parameters.For each month j and each mission M , the time series is constructed as These time series refer to the center of the cell and are corrected for backscatter related penetration effects.A schematic illustration of the results of this step is given in Fig. 4a.Otherwise, the offset with respect to the reference mission (Envisat) is obtained from the differences in the surface elevation parameter a0 of the fit according to Eq. (2) (dashed dark red lines; Seasat and Geosat).This is only applicable if all missions show a consistent rate of dh/dt as it is shown here.c) The specific offset between PLRA, SARIn and laser data depends on the sampling of the topography within each single cell.Nevertheless, the techniques can be combined by reducing each elevation time series by the specific elevation at the reference epoch (09/2010).To reduce the influence of non-linear surface elevation changes, this reference elevation is obtained from a 8-year interval only (gray area).

Merging different missions and techniques
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;Khvorostovsky, 2012).When merging data from different observation techniques (PLRA, SARIn 5 and laser) the calibration gets even more challenging.
We chose an approach in different steps which is depicted in Fig. 4. In a first step, we merge all the PLRA data.For these missions the topographic sampling of the instruments is similar and thus the offsets are valid over larger regions.For in Fig. 4b).Outliers in the regional varying offsets were removed using a moving median filter and the result is smoothed using a gaussian filter (σ = 20 km).Maps of the offsets with respect to Envisat are shown in the supplementary material Fig. S3.For the ERS missions, we find significant differences in the offsets for ice and ocean mode, hence, we determine separate offsets for each mode.The ERS-2 (ice mode)-Envisat offset distribution looks very similar to the pattern presented by Frappart et al.
(2016) 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, which do not have an overlap with other missions, we used the parameter a 0,M from Eq. ( 2).
In Fig. 4b the joint trend estimation is depicted by the equal slope of the dashed red lines.The different reference elevations a 0,M at t 0 = 09/2010 relate to the calibration offsets.This method, however, is only valid if the rate of elevation change is sufficiently stable.To find regions where we can assume such a stable rate, we compared the multi-mission dh/dt from Eq.
(2) with those, obtained from similar single-mission fits (Eq. 1) for ERS-1, ERS-2, Envisat and CryoSat-2 respectively.Only where all the differences between the single-mission rates and the combined rate are smaller than 5 cm/yr, we consider the disturbances due to interannual variations or an acceleration of the rate as not significant.Maps of these offsets with respect to Envisat are shown in the supplementary material Fig. S4.Regions where this stability criterion is fulfilled are mainly found on the plateau.The mean values over all cells amount to -0.85 m for Seasat and -0.72 m for Geosat.The corresponding standard deviations of 0.51 m and 0.34 m respectively are mainly a result of the regional pattern of the offsets.However, here we do not apply a regionally varying offset as we have no information how the offsets could be extrapolated towards the coastal regions where the stability criterion is not fulfilled.Furthermore, parts of the pattern might be explained by SEC anomalies prior to 1992 while afterwards the rate was sufficiently stable.Hence, applying the mean bias is the most robust solution.Within its error bars, it agrees very well with the ocean/sea ice based bias of -0.77 m between Seasat and Envisat (Fricker and Padman, 2012).With the help of these biases, all PLRA missions were corrected towards the chosen reference mission Envisat and thus form a joint PLRA time series (red in Fig. 4c with additional CryoSat-2 LRM mode where available).
To determine the biases of CryoSat's SARIn mode and ICESat compared to PRLA, simultaneously observed epochs could be used again.However, as noted in Sect.coherent results can also be obtained from the first 15 years by the early missions in these regions.From panels d and e, compared to 1992-2017, we see that the large persistent change rates, which are mainly related to ice dynamics, can also be observed over shorter time intervals.However we also see significant differences between the two intervals, leading to the conclusion that interannual variations or changing ice dynamics may have a strong influence on these time intervals.Hence, to study the temporal variability of the observed elevations, in the following we analyzed the full time series over different spatial and temporal scales.
At scales of subregions of the Antarctic Ice Sheet or drainage basins we calculated cumulative volume time series, shown in Fig. 7 and 8. Therefore, we multiplied the elevation time series by the corresponding grid cell area and summed them up over the respective region.We used a modified version of the drainage basin definitions by Rignot et al. (2011).This set was modified for an update of the ice sheet mass balance inter-comparison exercise (IMBIE, Shepherd et al., 2012) and is depicted in Fig. 6a.To correct for elevation changes due to GIA, all time series were corrected using the model IJ05_R2 (Ivins et al., 2013).This GIA-model predicts an uplift of 5 mm/yr near the Antarctic Peninsula and rates between -0.anomalies from ECMWF ERA-Interim (Dee et al., 2011) and GIA-corrected GRACE mass balances (Groh and Horwath, 2016) to our data.A direct comparison would impose the need to convert our volume changes to mass.This conversion, however, is 5 subject of controversial discussions whether a density mask (e.g.Sørensen et al., 2011;McMillan et al., 2014) is appropriate or if a firn densification model (e.g.Zwally et al., 2015;Kallenberg et al., 2017) is able to correctly model the interannual to decadal variations, especially in East Antarctica.As this is beyond the scope of this paper, all mass time series data are plotted at a scale of 1/3 of the volume scale.Hence, for a mass variation of fresh snow the curves would be directly comparable while for a mass loss or gain at a higher density, the mass curves vary stronger than the volume curve.14 The Cryosphere Discuss., https://doi.org/10.The spatial pattern of the elevation changes from year to year are shown in Fig. 9. Yearly averages are calculated from July until June of the next year and subtracted from each other.For 1978, the average is calculated from July to October only, the mission lifetime of Seasat.For the last period, denoted as "17", the average spans from July to December.The time intervals marked in yellow refer to differences that span more than one year.A similar plot for the precipitation anomaly data can be found in the supplementary material (Fig. S8).
To summarize the results we fitted trends to the volume change time series over Antarctica and several subregions which are listed in Tab. 2. It should be noted that, due to the different orbit inclinations of the missions, only the respectively covered regions can be used to calculate trends.The omissions of the ERS/Envisat polar gap can be assessed, at least for the period 2010-2017, by comparing the volume trends that respectively include and exclude the area south of 81.5°S.The difference (-151±14 km 3 /yr versus -170±11 km 3 /yr) is as low as 19 km 3 /yr.The error measure given with the rates here is the standard deviation from the linear fit.In the stochastic model, also a covariance between consecutive epochs has been considered but systematic sources of error could not be accounted for, hence our error estimate tends to be too optimistic.
Comparing the rates over different time periods reveals that fitting a linear trend is highly dependent on the time interval used.For the whole ice sheet north of 81.5°S we obtain an average volume rate of -34±5 km 3 /yr for 1992-2017.Separating the time series at 2010 when several glaciers started to accelerate, we obtain a slightly positive rate (27±5 km 3 /yr) for 1992-2010 but a dramatically larger volume loss of -170±11 km 3 /yr for 2010-2017.This again points out that linear rates are only comparable when they refer to the same time interval.
To validate our results, we used a set of 19 kinematic GNSS profiles (Schröder et al., 2017).They have been observed profiles was assessed by a crossover comparison between independent profiles of the same season.It is in the range of 4 to 9 cm.
One profile (K08C) has not been used as the poorly determined antenna height offset might impose larger errors.In Schröder et al. (2017) we used these profiles to validate the absolute altimetric elevation measurements.Here, we only work with elevation differences w.r.t. a reference epoch.As we want to analyze the temporal variability, we compare elevation changes over the same time observed by both techniques.For each crossover difference between kinematic profiles from different years we compare the differences of the respective altimetric SEC epochs in this location (δ∆h = ∆h KIN − ∆h ALT ). Figure 10 shows the results of this validation.We obtain an overall agreement of 6±16.5 cm or an RMS of 17.7 cm.As both techniques should observe the same elevation change, the difference δ∆h can be interpreted as an error measure.It can be expressed as as two kinematic profiles (with the respective RM S KIN ) and two altimetric epochs (with RM S ALT ) contribute to this difference.Using the a priori RM S KIN (4-9 cm) and resolving for RM S ALT yields the empirical RMS of the altimetric SEC data which lies in the range between 8 and 12 cm.A comparison with the satellite altimetry measurement precision (Tab. 1) shows that this corresponds to what we expect from the underlying measurements.The majority of the profiles was observed between 2007 and 2015.Hence, the differences also include different missions and observation techniques.So, we can conclude that no significant additional uncertainty has been added due to the combination.Small positive or negative biases in some regions of the map could be attributed to remaining antenna height offset errors in one of the GNSS profiles, but might also originate from remaining penetration effects of the radar signals.The very high agreement between the two datasets shows that our processing successfully eliminated the biases between the different satellite altimetry missions and thus provides reliable results.Nevertheless, only altimetric SEC within the interval 2001-2015 can be validated in this way.For the earlier missions, no spatially extensive high precision in situ data are available to us.The linear elevation change rates in Fig. 6 show the regions which experience a significant thinning (Amundsen Sea Embayment, Totten Glacier) or thickening (Kamb Ice Stream) which was already shown by a range of previous publications (e.g.Wingham et al., 2006;Flament and Rémy, 2012;Helm et al., 2014).In contrast to those studies our rates are calculated from the combined time series and hence, cover a much longer time interval of 40 years for the coastal regions of East Antarctica (Fig. 6a) and 25 years also for large parts of the interior of the ice sheet (Fig. 6b).In both of these long-term trends we see elevation gains of several cm/yr along the Antarctic Peninsula and in most of the coastal regions of East Antarctica between 45°W and 110°E.In Wilkes Land (C'-D') we see persistent rates of volume losses not only at Totten Glacier but also on some smaller glaciers as Frost, Mertz, Ninnis and Cook Glacier.Rignot (2006) observed previously that these glaciers, which are grounded below sea level, are thinning and losing mass.Our 25 year interval furthermore reveals the stability of the interior of the East Antarctic Ice Sheet.For 54% of the observed area of East Antarctica (north of 81.5°S), we find that the surface elevation changes are less than ±1 cm/yr.
The different sub-intervals in Fig. 6c-e show that the spatial pattern of elevation changes in East Antarctica is not constant over time.During 1978During -1992 the East Antarctic sector west of the Amery Ice Shelf (A-B) was losing volume while in the eastern sector (C-D') the rates were mainly positive.For the interval 1992-2010 we see positive rates in Coats Land and western Dronning Maud Land (J"-A') and in Princess Elizabeth Land (C-C') while elsewhere, the rates are very close to zero (except for the dynamically thinning glaciers).For 2010-2017 we see strong elevation gains in Dronning Maud Land and Enderby Land (A-B) but also strong losses in Princess Elizabeth Land (C-C').In the supplementary material Fig. S5 we calculated similar rates for the precipitation anomalies.These precipitation induced rates show a very consistent pattern compared to the elevation changes in East Antarctica.For the interval 2010-2017 they show that the strong mass losses in Princess Elizabeth Land can be related to a similar decline in precipitation.Also in Dronning Maud Land or Wilkes Land, the rates over the different intervals agree very well.In contrast, for West Antarctica, the patterns are significantly different which supports that the elevation changes are driven by ice dynamics here.Here, we can observe how the dynamic thinning has spread further inland during the last decade as described by Konrad et al. (2016).The strong variation between the different intervals points out that a linear rate over a decade may significantly differ from the long-term rate.
To cope with this interannual variation, we thus analyze the full time series of SEC, obtained after combining the singlemission time series as described in Sect.3.2.Figure 5a Wouters et al., 2015) and in the Getz and Abbot region (F-G, Chuter et al., 2017) started to loose mass at steadily growing rates.The individual basin time series which can be found in the supplementary material confirm this increasing loss.Since the study of Wouters et al. (2015), however, we observe that the region H'-I as well as the northernmost tip of the Peninsula (I-I") regained volume.Comparing the time series to anomalies in precipitation reveals that this can be explained by an extreme snowfall event in this area in 2016.The shape and orientation of the Peninsula makes GRACE observations challenging with respect to leakage and GRACE error effects (Horwath and Dietrich, 2009).Nevertheless, the results of the satellite gravity mission confirm this anomaly as well.Here, the maps of differences in yearly averages help to discriminate the spatial pattern of the anomalies.In Fig. 9 it becomes evident that an extreme snowfall event in 2016 occurred in western Palmer Land and northern Ellsworth Land, leading to this regain of volume.For the majority of the outlet glaciers in West Antarctica, the yearly differences in volume and precipitation do not agree very well.Instead, here, very persistent patterns of volume change can be observed, e.g. in the Amundsen Sea Embayment or at Kamb Ice Stream, which support the interpretation of a dynamic origin of these changes (Joughin et al., 1999;Pritchard et al., 2009).Furthermore, in Fig. 7b and Fig. S6 b-f the mass time series from GRACE change at significantly larger rates compared to the volume time series.This also indicates that changes occur at remarkably higher densities.
In contrast to West Antarctica, the time series summed up over the East Antarctic Ice Sheet (EAIS) is even slightly positive (see Fig. 7b).As it can be seen from the individual basins (Fig. 8  In this basin a difficulty in bridging the observational gap in the calibration of the old missions becomes evident.The altimetric elevations since 1992 show a linear behavior.Hence, we also use some data of this basin as described in Sect.3.2 to determine the mean calibration bias of Seasat and Geosat.The precipitation anomalies confirm the linearity since the early nineties but they also show a significant change in the trend before.However, by using a mean calibration bias over entire Antarctica, the impact of this presumable non-linearity in the mission calibration is likely compensated by many different biases in other regions. The differences of the yearly averages of SEC in Fig. 9 and the respective variations of precipitation in Fig. S8 reveal the spatial signature of the interannual variations.Even for the early missions Seasat and Geosat they show a very consistent picture and hence again demonstrate that these missions can provide important information to extend the observed time interval to a maximum.Between Seast and Geosat, we see elevations gains along the East Antarctic coast between 0°E and 70°E while west of the Amery Ice Shelf (70°E-150°E) the differences are negative.This is confirmed by the precipitation data for most of the area, even if the data here do not include 1978.Elevation gains in coastal western Enderby Land (40°E-50°E), which are very prominent in the precipitation data, are not well resolved by the early altimetry missions due to the distinct topography.
Nevertheless, we see some positive elevation changes there too.In 1986, we observe a positive anomaly in Princess Elizabeth Land and Wilkes Land which coincides with an increased snowfall.A very similar spatial pattern repeats again in 1996in , 2001in , in 2004in /2005in and in 2009in . Mémin et al. (2015) ) already observed such periodic elevation anomalies using Envisat and calculated a frequency of 4.7 years.They identify the Antarctic Circumpolar Wave as a main driver for this periodic signal but also note that it is superimposed by ENSO.This might be the reason why in 2013 another occurrence, which would by expected by the periodicity, is unusually weak.
In contrast to the highly variable coastal regions, the interior of the EAIS is very stable.As it can be seen from Fig. 7e the surface in the LPZ varied between 5 and -8 cm during the last 25 years.A very slightly positive rate is consistent with the precipitation anomalies over this period.Zwally et al. (2015) observed a positive rate in this region as well and explain it by a continuing dynamic response to increased accumulation since the early Holocene.However, their estimated rate of about 2 cm/yr would result in an elevation increase of 0.5 m over the past 25 years, which clearly contradicts our results.
Furthermore, kinematic GPS measurements on snowmobiles around Vostok station (Richter et al., 2014) yield an elevation change rate of 0.1 ±0.5 cm/yr for the period 2001 to 2013 which is in every good agreement with our altimetric time series.

Conclusions
In this paper we presented an approach to combine different satellite altimetry missions, observation modes and techniques.
The reprocessing of the conventional pulse limited radar altimetry guarantees that two fundamental steps in processing of radar ice altimetry, the waveform retracking and the slope correction, are handled consistently.Furthermore, we showed that the advanced methods used in this processing improved the precision by up to 50%, compared to the widely used standard products. 20 The Cryosphere Discuss.Besides the linear rates we were able to create monthly multi-mission time series which allow to study the evolution of the ice sheet volume over the full time span covered by any of the measurements.Peaks and steps in the altimetric SEC time series agree very well with the corresponding cumulated precipitation anomalies from ERA-Interim as well as the mass balance time series from GRACE.Note that we do not expect a perfect physical correlation between the three records.Making them fully comparable would require a sophisticated separation between surface and ice dynamic processes and the consideration of firn compaction (Ligtenberg et al., 2011;Li and Zwally, 2011).Both are beyond the scope of this paper.Nevertheless, the comparison of the different datasets shows, that our methodology successfully eliminated the biases between the different missions and observation techniques.We do not see any jumps in the time series, which could be interpreted as a calibration bias.
We conclude that this paper shows how to produce an altimetric elevation time series that is free of obvious artifacts of processing or calibration biases.This is essential to analyze the time series over the full time span of up to 40 years.

Figure 1 .
Figure 1.Spatial and temporal coverage of the satellite altimetry missions used in this study.The colors of the time bars denote the maximum southern extent of the measurements and thus the size of the respective polar gap.

TheFigure 4 .
Figure 4. Schematic elevation time series as obtained from Eq. (3) to illustrate the combination of the different data sets.a) Single-mission time series with inter-mission offsets.To emphasize the different techniques, pulse-limited radar altimetry (PLRA) observations are colored in red and pink, data from SARIn mode of CryoSat-2 in blue and the laser altimetry measurements of ICESat in cyan.b) To combine the PLRA data, overlapping epochs are used where they exist (ERS1-ERS2, ERS2-Envisat; red area).Otherwise, the offset with respect to the overlapping missions (ERS-1, ERS-2, Envisat, CryoSat-2 LRM) the offsets are calculated from simultaneous epochs (red area 8 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-49Manuscript under review for journal The Cryosphere Discussion started: 19 March 2018 c Author(s) 2018.CC BY 4.0 License.
3.1 and Fig.3, the different sampling of topography by the three techniques might lead to completely different surfaces, fitted to the respective elevation measurements.In Eq. (2), this is accounted for by technique specific slope parameters but also the reference elevation (a 0,M ) contains a component of the difference of the effective topography.Hence, the biases between different techniques refer to the effective topographies fitted to the measurements only and are distinct for every grid cell.Only cells which contain measurements from different techniques would hence allow to estimate the bias for this cell.Due to different orbit inclinations, this applies only to very few cells and consequently, the use of individual biases from simultaneous epochs is not satisfactory.Instead, we applied another linear fit to the three, technique specific SEC time series from SARIn, laser and the merged PLRA time series.We estimate a joint dh/dt and an individual reference elevation a 0,T for each technique T for the epoch 09/2010.As unmodelled effects like interannual variations or accelerations in the rate might adulterate the parameters, we restrict this fit to a short time interval of 8 years (09/2006-09/2014, gray area in Fig. 4c) which contains all techniques.The The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-49Manuscript under review for journal The Cryosphere Discussion started: 19 March 2018 c Author(s) 2018.CC BY 4.0 License.technique-specificreference elevations a 0,T are subtracted from the respective time series and hence, all series of elevation differences refer to 09/2010 and can be combined.This final combination of the techniques is performed using a weighted spatio-temporal averaging with 20 km σ gaussian weights in spatial domain and including the two consecutive epochs in the temporal domain to reduce the noise of a single pass.This averaging is used to interpolate unobserved cells as well.Due to the smoothing of the weighting function, we reduce our spatial grid resolution to 10x10 km now and did not interpolate a value to cells which are mainly covered by rocks(Burton-Johnson et al., 2016).To avoid extensive extrapolation to unobserved regions, the data around each cells center was classified into six sectors and interpolation was only performed if at least three of the sectors contained data.4ResultsWe obtain grids of surface elevation change (SEC) with respect to 09/2010 for each month observed and at a 10 km spatial resolution from the approach described above.Five examples given in Fig.5ademonstrate the spatial coverage and the ability to resolve the temporal evolution.The time series along the profile from Totten Glacier to Lake Vostok (Fig.5b) shows the consistency of the data, including the two early missions.Especially in the coastal areas where the signals are large, the early missions can provide important information to extend the time interval by additional 15 years.Surface elevation change rates fitted to the combined SEC time series over different time intervals are given in Fig.6.For the interval 1978-2017 these rates show which processes persist over four decades.Due to the orbit inclination of the early missions, these results are mainly limited to coastal areas in East Antarctica but also some flat regions along the Antarctic Peninsula are covered by observations of the early missions.A spatially more comprehensive picture is given in Fig.6b, covering the last 25 years since the launch of ERS-1.Similar rates calculated over sub-intervals of the full time range, namely07/1978-12/1992, 04/1992-12/2010 and 01/2010-12/2017, are displayed in the panels c-e.The interval 1978The interval  -1992   shows that

Figure 5
Figure 5. a) Five example snapshots of the combined surface elevation time series.The height differences refer to our reference epoch 09/2010.b) Yearly means of surface elevation change along the profile marked by the black line in a) from Totten Glacier towards Lake Vostok.

Figure 6 .Figure 7 .Figure 8 .
Figure 6.Multi-mission surface elevation change rates from the combined SEC time series over different time intervals.The solid lines mark the drainage basin outlines, the dashed line shows the outline of the low precipitation zone.a) Trends for the area covered by observations since 1978.b) Trends since 1992 for the area covered by the ERS missions and Envisat.c-e) consecutive time intervals show coherent patterns in dynamic regions but also large variations due to interannual variability.

Figure 9 .
Figure 9. Elevation difference between the yearly mean SEC (July to June of next year) of consecutive years.Marked in yellow are differences spanning more than one year due to altimetry data gaps.

Figure 10 .
Figure 10.Map (a) and histogram (b) of elevation differences observed at crossover locations of GNSS profiles minus the elevation difference obtained from the altimetric results for the respective times.
shows five months as examples of the spatial coverage for the respective epochs.The profiles in Fig.5bat different epochs show that in the southernmost part, which is located in the region of the subglacial Lake Vostok, the elevations vary by less than ±10 cm over the last 25 years.Around kilometer 600 where the profile bends into the main flowline of Totten Glacier, we see a significantly rising elevation.The profiles at different epochs reveal that this is not a continuous change but that there is a distinct jump in the early 2000s.From the yearly precipitation anomaly maps (Fig.S8) we can see that in 2001 a strong snowfall event occurred in this region.At the very front of Totten Glacier, the elevation dropped by almost 12 m from 1985 to 2010.Seasat could not successfully track a return signal near the front but the profile of 1978 shows that also before 1985, the glacier was already thinning.Between 2010 and 2017, the elevation decreased The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-49Manuscript under review for journal The Cryosphere Discussion started: 19 March 2018 c Author(s) 2018.CC BY 4.0 License.by more than 7 m.This means that the thinning accelerated from about 0.5 m/yr to more than 1 m/yr in the last decade.These rates are consistent with the approximately 0.7 m/yr Flament and Rémy (2012) obtained from Envisat (2002-2010).Li et al.(2015) observe a grounding line retreat at Totten Glacier of 1 to 3 km between 1996 and 2013 using SAR Interferometry.They conclude that this indicates an ice thinning of 12 m (or 0.7 m/yr) which is also very consistent with our results.Summed up over larger regions the measurement noise is reduced and the time series reveal the losses and gains over the respective area.As can be seen from Fig.7athe sum of the Antarctic ice volume change is slightly negative over the last 25 years.This, however, is mainly a result of the accelerated dynamic thinning, especially in West Antarctica, over the last decade.Konrad et al. (2016) reported a spread as well as an amplification of the thinning of the glaciers along the Amundsen Sea Embayment (basin G-H).Since around 2007 to 2010, also the adjacent glaciers along the Bellingshausen Sea (H-H', and S7), the losses over the last decade at Totten Glacier (C'-D) and in George V Land (D-D') are compensated by volume gains in the remaining basins (A-C').Here, the two early missions Seasat and Geosat help to extend the time series and thus get a better understanding of what is the long term behavior.Peaks and steps in the altimetry time series can be nicely related to corresponding events in the precipitation.We observe accumulation events in Dronning Maud Land (A-A') in 2002, 2009 and 2011.The yearly difference maps for SEC (Fig.9) and precipitation (Fig.S8) agree very well and show that, in 2002 the extreme snowfall occurred mainly in western Dronning Maud Land while in 2009 and 2011 the entire basin as well as the consecutive Enderby Land (A'-B) were affected.In Princess Elizabeth Land (C-C') and the catchment of Lambert Glacier, we observe the effect of a decreased precipitation in 1993 and 1994 which lead to a significant surface lowering.Conversely, in 2001/2002 a very strong snowfall event compensated these losses.Due to the large extend of these basins, this dent can even be detected in the time series over the entire Ice Sheet.Also in George V Land (D-D') the strong decrease in elevation between 2008 and 2010 is highly correlated to a deficit in accumulation.The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-49Manuscript under review for journal The Cryosphere Discussion started: 19 March 2018 c Author(s) 2018.CC BY 4.0 License.
, https://doi.org/10.5194/tc-2018-49Manuscript under review for journal The Cryosphere Discussion started: 19 March 2018 c Author(s) 2018.CC BY 4.0 License.From the combined time series of SEC we calculated the long-term surface elevation change rates over the last 25 years.In the coastal regions of East Antarctica we extended the time series back to 1978 using the observations of the missions Seasat and Geosat.Hence, we were able to calculate the elevation change rates covering four decades there.Comparing rates calculated over different time intervals revealed which signals are persistent and where the interannual variations are dominant.The map of long-term trends, obtained from an unprecedented time interval, shows that large parts of the East Antarctic plateau are very close to equilibrium.
Author contributions.L. Schröder designed the study and developed the PLRA reprocessing, the repeat altimetry processing and the time series generation.V. Helm supplied the reprocessed CryoSat-2 SARIn data.All authors discussed the results and contributed to the writing and editing of the manuscript.The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-49Manuscript under review for journal The Cryosphere Discussion started: 19 March 2018 c Author(s) 2018.CC BY 4.0 License.
Illustration 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, the planar surface fitted to the measurements (dashed lines) should not be mixed over different techniques.

Table 2 .
Volume change rate for Antarctica (ANT), the East and West Antarctic Ice Sheet (EAIS and WAIS) and the Antarctic Peninsula (APIS).The size of the observed area refers to observations from: Seasat for latitudes north of 72°S, ERS-1/2 for latitudes north of 81.5°S and CryoSat-2 for regions without a southern limit.The very sparse observations of Seasat and Geosat on the Antarctic Peninsula do not allow to calculate a reliable long-term trend.