Validation of satellite altimetry by kinematic GNSS in central East Antarctica

Ice-surface elevation profiles of more than 30,000 km in total length are derived from kinematic GNSS (GPS and the Russian GLONASS) observations on sledge convoy vehicles along traverses between Vostok station and the East Antarctic coast. These profiles have accuracies between 4 and 9 cm. They are used to validate elevation datasets from both radar and laser satellite altimetry as well as four digital elevation models. A crossover analysis with three different processing versions of Envisat radar altimetry elevation profiles yields a clear preference for the relocation method over the direct method of slope 5 correction and for threshold retrackers over functional fit algorithms. The validation of CryoSat-2 low-resolution mode and SARIn mode datasets documents the progress made from baseline B to C elevation products. ICESat laser altimetry data are demonstrated to be accurate to a few decimetres over a wide range of surface slopes. A crossover adjustment in the region of subglacial Lake Vostok combining ICESat elevation data with our GNSS profiles yields a new set of ICESat laser campaign biases and provides new, independent evidence for the stability of the ice-surface elevation above the lake. The evaluation of 10 the digital elevation models reveals the benefits of combining laser and radar altimetry.


Introduction
Surface elevation data are crucial for a broad range of applications in polar sciences.Only satellite altimetry is able to provide this information with a high and nearly uniform accuracy and precision for almost the entire Antarctic ice sheet.This high accuracy also allows us to infer temporal changes in ice surface elevation, which is of prime scientific interest in the context of ongoing climate change (Shepherd et al., 2012;Groh et al., 2014).However, systematic effects (as varying surface properties or measurement biases) can deteriorate the derived elevation trend ḣ and -if not corrected thoroughly -they might lead to misinterpretation of the observations (Arthern et al., 2001;Lacroix et al., 2009).
One crucial step in the processing of surface elevations from satellite radar altimetry (SRA) over ice sheets is the slope correction (Brenner et al., 1983).Due to the size of the beam-limited footprint of about 20 km in diameter the first reflection can originate from a location up to several kilometres away from the nadir point in a sloping surface.Different approaches exist to correct for this effect (e.g.Bamber, 1994;Roemer et al., 2007) but, as the corrections can exceed 100 m (Brenner et al., 2007), remaining model errors may introduce height errors of up to several metres.This is the major factor limiting the application of SRA in the steep and rugged coastal areas (Flament and Rémy, 2012).
Another issue when deriving ice-surface elevations from SRA data is the penetration of the microwave signal into the upper firn layers.This results in a mixed return signal consisting of surface reflection and volume reflection (Ridley and Partington, 1988).Here, the selection of an appropriate retracking algorithm is essential.One approach to minimize the influence of the volume echo on the observed surface elevations is to retrack at the very beginning of the waveform (Davis, 1997).Another method is to apply appropriate corrections using parameters of the radar waveform shape (Wingham et al., 1998;Flament and Rémy, 2012;Zwally et al., 2015).
For the Ice Cloud and land Elevation Satellite (ICESat) mission the effects of topographic correction and signal penetration do not arise or are negligible as the on-board altimeter uses laser signals.Hence, significantly higher accuracies can be achieved.Nonetheless, those measurements are also not free of systematic errors.Pointing errors and orbital variations (Luthcke et al., 2005) or saturation effects (Scambos and Shuman, 2016) may cause laser campaign biases which induce spurious trends of up to 2 cm yr −1 (Hofton et al., 2013;Gunter et al., 2014) and introduce errors of more than 100 Gt yr −1 in mass balance estimates (Hofton et al., 2013).
In order to quantify the impact of these errors and to evaluate methods for their correction, independent elevation data of high precision and accuracy is crucial.Here, we make use of ice-surface elevation profiles in central East Antarctica comprising more than 30 000 km in total length.These profiles are observed with kinematic GNSS (Global Navigation Satellite System, which means the Global Positioning System (GPS) and the Russian Global Navigation Satellite System (GLONASS) in this case) carried out over more than one decade on sledge convoy vehicles along continental traverses.
2 Surface elevations from kinematic GNSS profiles

Kinematic GNSS observations
The Russian research station Vostok is located in the central part of East Antarctica (106.8 • E, 78.5 • S).It is the main base for a wide range of scientific fieldwork related to the subglacial Lake Vostok.Between 2001 and 2015 several kinematic GNSS profiles have been measured in the area of the lake as well as on the scientific traverses from Vostok Station to the East Antarctic coast.
Geodetic dual-frequency GNSS receivers with external antennas were used for kinematic profiling as well as on the reference stations.Two different types of profiles can be distinguished with respect to the vehicles onto which the GNSS antennas were mounted.The first type are observations performed on lightweight snowmobiles.With the help of such profiles Richter et al. (2014a) have shown that the surface elevation around Vostok Station has been stable over the last decade, confirming the results of permanent GNSS observations (Richter et al., 2008(Richter et al., , 2014a)).Such profiles acquired on snowmobiles provide accuracies of only a few centimetres (see also King et al., 2009;Siegfried et al., 2011) and are thus well suited for precise studies on local elevation and elevation changes.However, due to logistic reasons they usually cover only a very limited area and are therefore not considered here.
The second type of observations are carried out on heavy convoy vehicles.Those are tractors on tracks, designed to pull sledges with containers for accommodation and fuel tanks (Fig. 1).Hence, they are ideal platforms for such measurements over very long distances.This is a precondition for the validation of satellite altimetry on a larger scale, as it helps to minimize the influence of regional peculiarities, especially due to specific topographic conditions (Kohler et al., 2013).The disadvantage of such heavy platforms, compared to snowmobiles, is that they sink into the soft upper snow layers up to several decimetres.The amount of the vehicle's subsidence, and thus of the height of the antenna above the snow surface, varies locally.Therefore, this antenna height has to be measured as often as possible along the traverse.
In the austral summer 2001/2002, our first surface elevation profile was acquired during a seismic convoy of the Russian Antarctic Expedition (RAE) along a 150 km transect in the southern part of Lake Vostok.During this traverse over 6 days a GNSS antenna was installed on the roof of a trailer, pulled by a traverse vehicle.
Starting from 2006, significantly longer profiles have been measured.In that time Mirny Station (93.0 • E, 66.6 • S) was the coastal logistical hub for the supply of Vostok Station by overland traverses.Several scientific observations were performed during these convoys (Masolov et al., 2001;Richter et al., 2013;Popov, 2015;Ekaykin et al., 2017).In the austral summer 2006/2007 kinematic GNSS profiles between Vostok and Mirny were observed on two convoy vehicles.For the first time these profiles cover the entire distance of about 1600 km, from the remarkably flat ice surface above Lake Vostok down to the rugged terrain at the coast.In the following season these profiles were repeated on two vehicles.
In 2009, Progress Station (76.4 • E, 69.4 • S) became the main logistic hub for Vostok Station.Further, an initial reconnaissance traverse from Progress to Vostok and back was performed in 2007/2008 that already included geodetic GNSSequipment.Since austral summer 2009/2010, several profiles between these two stations were measured each season.A number of different routes were used according to the needs of the participating scientific groups, snow conditions or logistical constraints.Figure 1 gives an overview of the locations and timing of the routes as well as two examples of the types of vehicles used.

GNSS data processing
We used the Bernese GNSS Software 5.2 (Dach et al., 2015) for the differential post-processing of the kinematic observation data.In such remote locations, other studies (Siegfried et al., 2011;Kohler et al., 2013) used the precise point positioning (PPP) technique, which does not require reference stations.However, PPP depends on precise high-rate satellite clock information, which is not available for our early campaigns.Geng et al. (2010) showed that both techniques are able to reach very high accuracies (∼ 3 cm) for very long distances to reference stations.The processing provides a 3-D coordinate of the GNSS antenna for each observational epoch in the terrestrial reference frame IGS08.Using those coordinates, a profile of ellipsoidal elevations referenced to WGS-84 is derived.For many of the profiles multi-system GNSS receivers were used; i.e. observations from the Russian GLONASS are logged in addition to GPS.The increased amount of observations improves the reliability of the solution significantly.As kinematic reference sites in this differential positioning we utilized static observations from campaign sites, for example in Mirny or Progress, from an own permanent receiver in Vostok installed in early 2008 (for details see Richter et al., 2014a) and additionally from the sites Casey and Davis of the IGS-network (see Fig. 1).To cope with the scarce 30 s sampling interval of the IGS-sites, those static observations had to be interpolated to the rate of the kinematic receivers (mainly 5 or 15 s, sf.Table 1).
For that purpose we used WaSoft, a software tool developed by Wanninger (2000).We adopt the processing strategy of Fritsche et al. (2014), which corrects or parameterizes the tropospheric and ionospheric delay, the antenna phase centre offsets and variations, solid earth tides and loading displacements.Special attention is paid to the resolution of the GNSS carrier phase ambiguities of the differenced observations.When the vehicle is halfway between Vostok and the coast, no baseline to a static reference station is shorter than 800 km.Then, only very robust ambiguity resolution strategies are able to produce satisfactory results.Therefore, in this case we used the Melbourne-Wübbena and the quasiionosphere-free linear combination only (Dach et al., 2015).As this is the most critical step in processing, a thorough outlier screening of the fixed ambiguity solutions is essential.
For this reason, we always used more than one, but typically four to five, baselines to different reference sites.Undetected cycle slips lead to very large deviations of the affected baseline.Thus, by processing each baseline independently and comparing the results to the combined solution, the baseline causing large deviations can be identified and the undetected cycle slip has to be introduced manually.

Derivation of surface elevation profiles
The antenna trajectory resulting from the GNSS positioning has to be corrected for the height of the antenna above the local snow surface in order to derive surface elevation profiles.This vertical offset is not constant as the amount to which the vehicle sinks into the snow depends on the regionally varying surface snow properties, but also on the vehicle type (e.g.track width).For example, Kohler et al. (2013)  to employ additional laser measurements on another vehicle in order to retrieve the amount of vehicle subsidence because the antenna height was not measured repeatedly along their profiles.During our traverses we measured this antenna height offset AH several times for each observation day and for each profile.However, the representativity of a single offset measurement may still be limited due to small-scale surface structures (sastrugi) at the locations of the measurements.Thus, to obtain a specific offset for each single epoch i, we use a regional average where d −1 ij is the inverse distance between the position at epoch i and the position of the antenna height measurement AH j .The offset of profile K08C was measured only once.
Here we model the subsidence from similar profiles with comparable vehicles.
Furthermore, a permanent tilt of the moving vehicles had to be considered.While driving in soft snow, especially when pulling heavy sledges, the nose of the vehicle gets lifted up while the rear buries deeper.This dynamic effect is not determined directly as our offset measurements are taken during stops when the vehicle stands upright.However, an instantaneous jump in antenna elevation from GNSS positioning is observed whenever the vehicle stops.Depending on the antenna's position on the vehicle, it can reach 20 cm.These jumps are used to correct the measurement AH for the vehicle dynamics.For this purpose we interpolate the elevations in movement (i.e.velocity > 1 km h −1 ) to the position of the antenna height measurement by fitting a quadratic function within a distance of 100 m around this point and comparing it to the average elevation in rest.
To reduce the noise and the influence of sastrugi, we applied a low-pass filter to the original antenna trajectory.Depending on the sampling interval (Table 1) and the velocity of the vehicles (about 7 km h −1 ), the usual point distance is in the range of 10 to 30 m.We applied a Gaussian filter with a Gaussian sigma of 60 m and a total length of 180 m.In addition, the trajectory positions are thinned out to an equidistant interval of 30 m.This substantially reduces the data amount (e.g. when data were logged during overnight stops of the convoy) without loss of information.
A comparison between kinematic GNSS profiles and satellite altimetry products requires consistency in the reference system used (King et al., 2009).ICESat (Schutz and Urban, 2014) and CryoSat-2 (Schrama et al., 2010) refer to ITRF08.This is identical to IGS08 within sub-millimetre level (Rebischung et al., 2012).Envisat GDR-C orbits refer to ITRF05 (Cerri et al., 2011) but as this affects the differences in the order of some millimetres only (Rebischung et al., 2012), it is considered negligible here, too.However, the treatment of the permanent tide in the reference systems underlying both techniques has a significant influence here.According to McCarthy and Petit (2004) the International Terrestrial Reference Frame (ITRF) and consequently IGS08, is a conventional tide-free frame.Hence, all tidal effects including the permanent effect have been removed from the coordinates of the reference stations.Our elevation profiles are henceforth also conventional tide free.Altimetric elevations, in contrast, refer to the mean tide system.The GNSS elevations are converted to this mean tide system using Eq.(7.14a) of Petit and Luzum (2010), which is a function of latitude and amounts to about −10 cm at 70 • S.

Accuracy and precision
Precision estimates for the epochwise coordinates consist of estimates of the quality of the antenna positioning and an additional uncertainty due to the reduction to the snow surface.In a first step we assess the quality of the GNSS processing.The formal coordinate errors reported by the processing software are too optimistic as they do not account for non-white-noise components.A more realistic measure is found by comparing multiple baseline solutions.As mentioned in Sect.2.2, the ambiguity resolution is a critical step in the GNSS data processing.Unrecognized cycle slips can distort profile sections over several kilometres and are thus not removed by the low-pass filter.However, such instances are identified within independent solutions using different reference sites.The average baseline coordinate differences RMS BL are used to derive realistic estimates for the precision of the kinematic positioning.Mean baseline differences in the vertical component are shown for each profile in column RMS BL of Table 1 and are in general of the order of a few centimetres.
An additional source of uncertainty is imposed by the reduction of the GNSS antenna elevation to the snow surface.This reduction varies regionally due to varying snow surface characteristics.Thus, besides the error of the offset measurements themselves, the offset corrections, obtained by Eq. (1), contain an additional interpolation error.We assess both types of errors using semivariograms.Here, we fit a linear function to the squared differences between the measurements of the antenna height offsets, with respect to the distances between those measurements.We obtain a constant part of 6 cm, which relates to the uncertainty of the antenna height measurement itself.It is potentially affected by local surface features (sastrugi) and residuals in the dynamic tilt correction.The distance-related additional uncertainty is 0.25 cm km −1 and accounts for the specific distance between the location of the respective offset measurement and the location to be interpolated.Using these values, the precision of the antenna height reduction through inverse distance interpolation (RMS AH ) is derived.Hence, the total precision measure for a single surface elevation observation is obtained by (2)

is a combination of GPS and GLONASS
To account for possible errors in the subsidence modelling for profile K08C, we add additional 10 cm to the uncertainty there.
A rigorous empirical test for the absolute accuracy estimate is performed by the calculation of height differences at crossover locations of two different profiles of the same season.The elapsed time between the two passes over this location is typically between a few minutes and some days, thus the surface elevation is assumed unchanged.We consider those profiles to be highly independent measurements as the antenna/snow-surface offset is determined independently, the satellite constellations are usually completely different and the equipment used is different.During the processing the tropospheric correction has also been estimated individually for each profile.We assume the GNSS technique itself to be practically unbiased as it is used to define the IGS08 reference system as well.As the differences h are calculated from two passes, the accuracy RMS X of a single profile at the crossover location is given by For our profiles we obtain RMS X in the range of 4-9 cm (see Table 1).This is slightly higher than the RMS S because it also includes the effect of vehicle dynamics.Nevertheless, it is a conservative estimate as in these crossovers the eleva-tions of the second profile are affected by the disturbances of the snow surface originating from the first vehicle pass.
As an additional independent validation, we compare our kinematic GNSS profiles to airborne elevation measurements from Operation ICEBridge (Studinger, 2014).The region under investigation is very sparsely covered by flights, but on 26 November 2013 an Airborne Topographic Mapper (ATM) lidar profile crossed Lake Vostok twice.Brunt et al. (2017) compared ATM measurements to ground-based GPS profiles in Greenland and found biases between −11 and 7 cm with precisions of less than 9 cm.After applying the correction for the mean tidal system, we calculated crossover differences h ATM between the nadir measurements and our profiles (see Table 1, last column).The amount of crossover points of this validation is very low (maximum 9 with K07A), but still valuable to check RMS X as an estimate for the overall absolute accuracy.The average offset from the 39 crossovers with all profiles is 4.9 cm and the standard deviation of the differences is 10 cm.Hence, we find that our accuracy estimates are realistic.There are several other airborne profiles with a laser altimeter crossing the convoy route to Mirny, but with an average of 14 cm and a standard deviation of 36 cm, the crossover differences with these laser altimeter profiles are not adequate for such a comparison.
L. Schröder et al.: Validation of satellite altimetry by kinematic GNSS in central East Antarctica

Elevation changes
The observation times of our kinematic GNSS profiles do not exactly coincide with the periods of the satellite missions validated here.Therefore, when comparing the GNSS-derived elevations with altimetry products, it is crucial to know to what extent elevation changes occurred between their respective observation epochs.While crossover differences within one expedition are used for accuracy estimates, the elevation differences in crossovers between profiles of different years allow us to assess temporal rates of surface elevation changes ( ḣ). Figure 2a shows that the obtained surface elevation rates are very small over the whole area.In Fig. 2a they are averaged at 20 km blocks to reduce random noise.However, the rates shown may still be affected by systematic errors effective over longer distances.One potential error source, in addition to those mentioned in Sect.2.4, is the impact of human activities on the snow surface.The immediate vicinity of the stations is obviously heavily affected but this is not the only region which had to be handled with care.For 5 decades the convoy between Mirny and Vostok used the same route.Especially above 3000 m, the heavy convoy vehicles and cargo sledges had followed exactly the same track in order to cope with the soft snow.This resulted in enhanced snow compaction along the track and the accumulation of a continuous ridge of several decimetres in height which is even visible on satellite imagery.The elevations and elevation changes along this part of the traverse are not representative of this area and are thus excluded from all subsequent studies.The profiles acquired on the modern Kässbohrer tractors are not prone to this effect since their wider tracks and relatively small weight relieves these vehicles of the need to reuse a pre-existing track.
The largest rates, but also the largest variations, are found on the traverse to Mirny.In the lower-elevation parts this is not an effect of anthropogenic disturbance.As the snow is much harder there, the tractors do not repeat the exact tracks of their predecessors.Nevertheless, the rates obtained in this area rely solely on 1 year time spans between the measurements of those profiles and must therefore be treated with care.In the areas where the time span is longer, very small rates are obtained.Averaging all 18 000 ḣ by introducing weights according to the time span results in a mean elevation change over the entire area of 4 cm yr −1 and a standard deviation (σ ) of a single crossover rate of ±11 cm yr −1 .
A detailed look into the results in the Lake Vostok region is given in Fig. 2c and d.In order to avoid the limitations arising from short observation time spans, we used only crossovers spanning 5 years or more.The average h of the 492 crossover differences in this area is −0.1 cm yr −1 with a standard deviation of a single rate of ±2.4 cm yr −1 .The ḣ values are not uncorrelated, especially due to possible systematic biases (e.g.antenna height reduction) which affect multiple crossovers of a profile.Therefore, for the accuracy of the mean h we only consider the number of combinations of independent profiles (27) in the estimation, resulting in a standard error of ±0.5 cm yr −1 .This comprises both real variations in surface elevation change rate and observational uncertainties.These results agree very well with the elevation changes observed by measurements on snowmobiles around Vostok Station (Richter et al., 2014a, 0.1 ± 0.5 cm yr −1 ).The latter profiles have a higher accuracy but yield a smaller amount of crossovers which may be affected by higher spatial correlation.
We conclude that the elevation change rates are very small in the region under investigation, but their spatial pattern is not determined with homogeneous reliability due to the short observation time span in some areas.For the following validations we consider the elevation change to be negligible.On the one hand we have chosen only missions which overlap in time with our profiles and thus the elevation changes due to the time differences are fairly small.On the other hand, at coastal regions where the rates might be larger, the errors of SRA are significantly larger too (metres for CryoSat-2 in SARIn mode, tens of metres for Envisat; see Sect.3.3).As this is not the case for ICESat, elevation changes will be discussed further in Sect.3.3.3.
3 Validation of satellite altimetry 3.1 Data

Envisat
We validate different altimeter missions to reveal characteristic effects of their respective techniques and approaches to derive optimum results.Conventional altimeter systems usually use signals in Ku-band and have a beam limited footprint of 10 to 20 km.Over ice sheets their signals penetrate into the upper firn layers.As an example for a conventional pulse-limited radar we validate the Envisat mission, operated by the European Space Agency (ESA).We use the Ku-band measurements of its altimeter system RA-2 acquired during the entire operation period (May 2002 to April 2012).In the Level 2 product (SGDR V2.1) the slope-induced error is corrected for using the relocation method (Bamber, 1994).This algorithm is designed to locate the measurement to the position where the first return signal comes from.The ESA data set contains results from two types of radar waveform retrackers applicable over ice sheets, ICE1 (based on the offset centre of gravity (OCOG) retracker by Wingham et al., 1986, with a threshold of 30 %) and ICE2 (a functional fit developed by Legrésy et al., 2005).We use both retrackers and compare their performance.
The Goddard Space Flight Center (GSFC) developed an in-house processing chain for radar altimetry with slightly different approaches for deriving surface elevations.Here, the direct method of slope correction was applied, which corrects the measurement at the nadir position.This reprocessed Level 2 data set is called the Ice Data Record (IDR) and also contains different retrackers.We use the GSFC V4 β-retracker as Brenner et al. (2007) summarize that this algorithm provides more accurate absolute elevations than threshold-based methods.
To remove potentially corrupted observations from the data, we used the measurement confidence flags (which are identical in the ESA SGDR and GSFCs IDR data sets) to find recorded distances out of range and to identify problems in the onboard processing and data handling, the ultra stable oscillator, the automatic gain control (AGC) or in the waveform samples.In addition to these instrumental errors we removed shots where the GSFC retracking algorithm failed as indicated by the retracking problem flag.In the SGDR data we furthermore used the fault identifier and the flag indicating that the ICE1 retracking in Ku-band was not successful.

CryoSat-2
Compared to the conventional SRA, ESA's CryoSat-2 has an improved resolution and accuracy due to its innovative design.In the smooth interior of the ice sheets the altimeter operates in the low-resolution mode (LRM) which is a conventional pulse-limited observation mode as in the missions before.Above steeper terrain, the altimeter is switched to SARIn mode.In this mode, the synthetic aperture radar (SAR) processing considerably improves the along track res-olution utilizing the Doppler/delay shift.Hence, the beamlimited footprint is subdivided in flight direction into stripes of roughly 250 m in length.The interferometric processing of the reception times at the two antennas allows the determination of the across-track angle to the point of closest approach (Wingham et al., 2006a).
We compare two different processing versions, Baseline B and C, of ESA's L2I data set.The "I" in the product identifier stands for the in-depth data set.It provides more parameters and flags and, over land, offers an additional feature relevant for our study.In the basic L2 product the SARIn ambiguity flag indicates an elevation difference between altimetry and a DEM exceeding 50 m.In this case the interferometric angle is considered erroneous and the measurement position is set to nadir.In the L2I product, however, this is not applied.This product allows us, therefore, to also validate the data at the margins where the a priori DEM itself is prone to large uncertainties (see Sect. 4).As an alternative approach, to identify outliers in the interferometric angle, we used the coherence flag and additionally excluded all measurements with a across-track angle exceeding 1 • (corresponding to the very edge of the antenna beam).Furthermore, we exclude all data where the respective retracker height error flag indicates problems in the determination of the retracking point.For Baseline B, the waveform is processed using the CFI retracker (Wingham et al., 2006a).In Baseline C two additional retrackers have been applied to the LRM data: a thresholdbased OCOG-retracker and another functional fit retracker www.the-cryosphere.net/11/1111/2017/The Cryosphere, 11, 1111-1130, 2017 called UCL, which is based on the Brown-model (Brown, 1977).

ICESat
In contrast to radar altimeters, the ICESat laser altimeter mission has a ground footprint of only 65 m and does not penetrate into the snowpack.Hence, surface elevation accuracies at the decimetre level can be achieved (Fricker et al., 2005;Shuman et al., 2006) which are almost comparable to our kinematic GNSS profiles.We use GLA12 elevation data (Zwally et al., 2014) from Release 34 (R34).We apply the saturation correction (Fricker et al., 2005) to the elevations and exclude all data where flags indicate off-nadir operation, orbit manoeuvres or any other factors degrading the orbit accuracy.We also remove data where the attitude flag indicates any problem with star trackers, gyro or the laser reference sensor.In order to exclude data affected by forward scattering in clouds or drifting snow (e.g.Siegfried et al., 2011), we reject all returns with a gain value exceeding 200, with a reflectivity below 10 %, with a misfit between the received waveform and the Gaussian model exceeding 0.03 V or where more than one waveform is detected (Bamber et al., 2009).

Crossover comparison
We validate the surface elevation data derived from satellite altimetry by applying the crossover method, outlined in Sect.2.4, to the intersections of the along-track altimetric profiles with our kinematic GNSS profiles.Therefore we interpolate the elevation of both profiles linearly if the distance between the two adjacent data points is less than 500 m.The along-track data point spacing amounts to 172 m for ICESat (Schutz et al., 2005), 300 m for CryoSat-2 (Wingham et al., 2006a) and 400 m for Envisat (ESA, 2007).In fact, the altimetry measurements represent average elevations over the respective footprint area.As we did not measure 2D-grids but straight lines, this cannot be taken into account here.The smoothing of the profiles with a filter length of 180 m, however, resembles the ICESat footprint at least in profile direction.
The largest error source for radar altimetry over a distinctive topography results from the slope correction.Brenner et al. (2007) found that crossover differences between ICESat and Envisat are less than 3 m for slopes below 0.1 • , but up to 50 m and more for slopes above 0. such a high sampling of slope.Instead, we investigate regions of different characteristic slope in separate histograms.This allows us not only to calculate a mean and standard deviation for each zone but also to identify deviations from a Gaussian distribution.Those histograms display the full range of results including potential outliers.In order to reduce their impact, an iterative 5σ filter is applied in the calculation of the mean and standard deviation.
We subdivide the region under investigation into four zones according to their mean surface slopes: > 0.5, 0.5-0.15,< 0.15 • and, as a subset of the latter characterized by extremely little surface roughness, the hydrostatic equilibrium area of subglacial Lake Vostok.The crossover differences between Envisat data and the GNSS profiles (Fig. 3) clearly demonstrate the relationship between surface slope and SRA errors and motivates our subdivision.The first zone comprises the coastal areas.Outliers and large errors in the SRA elevations are frequent there due to the rugged topography.On the other hand, it is the zone where the largest elevation changes would be expected.Hence, this zone introduces the largest uncertainties in ice-mass balance estimates based on SRA (Wingham et al., 2006b).The subsequent zone of intermediate slopes is still close to the coast and of low elevations.It may therefore also be subject to significant elevation changes.At the same time, SRA provides a higher accuracy there compared to the first zone.The third zone comprises the flat interior of the ice sheet.Here, elevation changes are generally small but, because of its vast areal extent, nevertheless important for mass balance studies.The ice above Lake Vostok, constituting the fourth zone, floats in hydrostatic equilibrium (Ewert et al., 2012).Surface gradients are very small and homogeneous in this area.Thus, the influence of the slope-induced error vanishes, offering a unique opportunity to study other effects such as the surface penetration of the radar signal.

ICESat campaign biases
Due to its smaller footprint size, the ICESat surface elevations are less sensitive to surface slope compared to SRA.However, the Geoscience Laser Altimeter System (GLAS) altimeter was operated in several laser operation campaigns due to laser degradation.Between the campaigns systematic biases exist (Fricker et al., 2005;Gunter et al., 2009).If not accounted for carefully, any systematic biases between campaigns can corrupt the inference of temporal surface elevation changes ḣ.To determine those biases, different surface types have been used, including the salt flat Salar de Uyuni (Fricker et al., 2005), the global oceans (Urban in Scambos and Shuman, 2016;Gunter et al., 2009), the ice surface above Lake Vostok (Ewert et al., 2012), the Antarctic lowprecipitation zone (Hofton et al., 2013;Gunter et al., 2014) or leads and polynyas in sea ice areas (Zwally et al., 2015).The estimated biases differ significantly between different data releases and, within the same release, depending on the surface type used for calibration.
Within the region around subglacial Lake Vostok (Fig. 1b, 100-108.5 • E, 76-79 • S) Ewert et al. ( 2012) applied a least squares adjustment of crossover differences between elevation profiles of ascending and descending ICESat (I ) orbits, acquired during laser campaigns i and j ( h I −I ij ).To cope with the lack of an absolute reference, these authors introduced a zero-sum condition.As a consequence, the laser campaign biases were determined as relative biases, relating to their overall average.This method relies essentially on the assumption of a stable surface throughout the ICESat observation period.This assumption is justified by the observational results of Richter et al. (2008).
Using the same region of ice-surface around Lake Vostok we derive a new set of laser campaign biases b for release 34.In addition to the ICESat crossovers (I ) used by Ewert et al. (2012) we also include crossover differences between ICE-Sat and our kinematic GNSS profiles ( h I −K iq ) and crossover differences between different GNSS profiles ( h K−K pq ).Including the unbiased GNSS profiles allows us to solve for the surface elevation change rate ḣ between the respective observation epochs t as an additional parameter and thus to overcome the assumption of a stable surface.We are hence able to separate real elevation trends ḣ from the apparent trend ḣb implied by the laser campaign biases.Furthermore, the incorporation of the unbiased GNSS elevations allows us to avoid the zero-sum condition and thus to determine absolute laser campaign biases.Combining all crossover differences results in three different types of observation equations: To account for the individual uncertainty of each GNSS profile, we introduce the epoch-wise surface elevation uncertainty RMS S as weights for the elevations.As shown in Sect.2.4 the empirical intra-expedition crossover differences RMS X are about 5 cm larger.For this reason we add 5 cm to each RMS S .For the ICESat elevations, we use campaign-specific average RMS X (Table 3) obtained from intra-campaign crossovers in the Lake Vostok region.

Envisat
The validation of the Envisat data (Fig. 4) shows that in the flat interior all processing versions provide precise elevations.Nevertheless, the crossovers over Lake Vostok reveal significant differences between the three retracker versions.
The standard deviations of the two functional fit retrackers yield similar results (51 cm for the β-retracker, 44 cm for ICE2).In contrast, the precision or the ICE1-retracked data is better by a factor of two (22 cm).This confirms the findings of Davis (1997), which argued for the superior precision of threshold retrackers.With respect to the kinematic GNSS profiles, the mean bias of all processing versions is negative.This can be explained by the penetration of the radar signal into the upper firn layers.However, significant differences between the retrackers are evident here.Compared to the ICE1-retracker, the mean reference surface of the ICE2 functional fit is 90 cm lower.Thus the influence of variations in firn pack properties is much stronger for this retracker.Between different data sets, a comparison of the biases of different retrackers should be treated with care.Here, elevation differences might also be caused by other factors as a different instrumental calibration value or alternative models for range correction.Compared to Lake Vostok, where the slope effect is negligible, significant differences can already be observed in the zone of least slopes (< 0.15 • ).Even the smooth topography there introduces additional uncertainties of about 30 cm for the relocated ESA data and 1.5 m for the GSFC data corrected by the direct method.Furthermore, the mean biases are shifted in the positive direction.The histograms reveal that this is a consequence of a deviation from the Gaussian distribution of the crossovers.The increased amount of positive differences means that the GNSS elevation is lower than the altimeter value.This is a consequence of the inability of the radar signal to observe depressions which are significantly smaller than the beam-limited footprint diameter (Brenner et al., 2007).For all versions, this effect increases progressively as the slope and hence the magnitudes of the L. Schröder et al.: Validation of satellite altimetry by kinematic GNSS in central East Antarctica depressions also get larger.In the intermediate surface slope zone the standard deviations of the ESA data sets grow to 3 m and in the coastal zone they grow up to 10 m.Here, the differences between the two ESA retrackers become negligible.For the GSFC product these errors are significantly larger.Hence, our results support a clear preference for the relocation method.Nevertheless, in zones of larger slopes the error of the slope correction becomes the dominating uncertainty contribution.
To avoid loss of tracking, Envisat switched the tracking bandwidth from high resolution to two lower-resolution modes when approaching steeper terrain.We analysed the performance of modes separately.Table 2 shows the results of each mode in the different zones.The number of crossovers indicates that the majority of data, even close to the coasts, was acquired in high-resolution mode.In the central zone the accuracy of the lower-resolution data is, as expected, worse.In the coastal zone, however, the other modes yield better results (if their small number of crossovers is considered as representative).The variations of the mean biases demonstrate, in turn, that these mode switches induce offsets in the data.We agree therefore with Brenner et al. (2007) not to use the sparse data of the lower-resolution modes for precise elevation change studies.

CryoSat-2
The results of the validation of CryoSat-2 elevations are shown in Fig. 5.A comparison with those from Envisat (Fig. 4) clearly shows the advantage of the SARIn mode in zones of larger slopes.Furthermore, a comparison of the recent Baseline C version with the previous Baseline B documents the improvements made by solving several issues.
A primary issue solved from Baseline B to C was a range bias of 67 cm in SARIn and 20 cm in LRM data (Scagliola and Fornari, 2015).In the SARIn data this bias reduction is clearly visible in all zones.Besides that, no major changes are evident and the standard deviations remain the same.Comparing the Baseline B LRM with the respective CFI-retracked version of Baseline C, we find a significant improvement in standard deviation.The refinements in the retracking procedure itself, described by Bouffard (2015), are probably responsible for those improvements.This holds true especially in the zones of stronger slopes where the retracking is more challenging.In the absence of slope-related effects on Lake Vostok, the correction of the range bias is also evident in the LRM data.
The main improvement in the performance of the Baseline C LRM data has been introduced by adding two additional retrackers (UCL and OCOG).The OCOG retracker shows standard deviations of about 20 cm over Lake Vostok, which is similar to the corresponding ICE1 retracker for Envisat.In contrast, the functional fit models show standard deviations of ∼ 50 cm which is similar to the results of the ICE2 retracker of Envisat.For the entire low-slope zone (< 0.15 • ) we obtain similar results when comparing CryoSat-2 LRM to Envisat.In the intermediate zone (0.15-0.5 • ) the two missions cannot be compared directly as the statistics for CryoSat-2 only relate to the subzone where the LRM is applied, which covers only the gently sloping areas.It should be noted that, even though SARIn mode is usually applied in coastal regions only, there are still some SARIn data available over Lake Vostok.On 28 July 2010 and the first week of June 2013 CryoSat-2 observed whole profiles across Antarctica in SARIn mode and also passed our region under investigation (including Lake Vostok) several times.These profiles allow us to directly compare the different modes and, in the case of Lake Vostok, the performance of their retrackers.The first column of Fig. 5 shows that the accuracy of SARIn is quite similar to the CFI retracker in LRM.
Different observational techniques have substantiated the stability of the surface elevation above Lake Vostok over timescales of typical satellite altimeter mission life times (Richter et al., 2008(Richter et al., , 2014a)).This stability, together with the low precipitation and the continuous monitoring of relevant parameters at Vostok Station, makes Lake Vostok an ideal area to examine apparent elevation variations in the altimetric time series.Spurious variations can be related to changes in surface backscatter and thus the backscattered power σ 0 of the altimetric signal (Wingham et al., 1998).Commonly, the relationship is determined as a regression coefficient and its influence is removed from the elevation time series (Wingham et al., 1998;Davis and Ferguson, 2004;Zwally et al., 2015).Figure 6 displays monthly averages of the crossover differences between the kinematic GNSS profiles and different CryoSat-2 LRM data sets.The Baseline B product (panel a) yields a high correlation as well as trends of the opposite sign in the elevations and the backscatter values.A trend derived from this elevation data set suggests a surface increase of 9.8 cm yr −1 which clearly contradicts all results from other studies.For the three retrackers of the Baseline C product (Fig. 6b-d) none of the backscatter curves show a significant trend any more.Nevertheless, the retracking methods based on functional fits (CFI in b, UCL in c) still exhibit a high correlation between h and σ 0 .Here, the retracking point is defined by the fit of the functional model to the whole waveform.Hence, it is more affected by volume scattering.In contrast the OCOG retracker uses a 25 % threshold of the OCOG amplitude and thus locates the retracking point much closer to the first radar return.The results in Fig. 6 indicate that threshold retrackers (panel d) produce the most precise elevations, especially in terms of repeatability.This confirms similar findings by Davis (1997).The seasonal variation of the signal almost vanishes.However, there is still a very small remaining amplitude, which correlates quite well with σ 0 .This indicates that there might still be some remaining effects of the snowpack properties superimposed on the elevation time series.Once the large variations disappeared, some jumps of a few decimetres are revealed.Apparent elevation jumps in two winters (2011/2012 and 2013/2014) Figure 4. Histograms, means and standard deviations of crossover differences between different Envisat data sets and kinematic GNSS profiles for four zones of characteristic surface slope.The range of the histograms is adjusted according to the values found in each zone, the colour scale is the same for all histograms.The displayed crossover differences contain uncertainties in the kinematic GNSS profiles (4-9 cm) and possible elevation changes between the observation epochs of both techniques in addition to the uncertainty in the Envisat data.
Table 2. Statistics of crossover differences between different Envisat resolution modes with kinematic GNSS profiles similar to Fig. 4. Outliers (> 5σ h ) are excluded iteratively.Each set of statistics contains h±σ h and the number of valid crossovers (sums can differ from the total number due to outlier rejection).Italic values are considered to be not statistically significant.to abrupt backscatter increase at the same time.Lacroix et al. (2009) detected a similar jump in Envisat data and referred it to changes in snowpack properties due to strong wind.However, the meteorological records from Vostok Station (wind, precipitation, temperature; not shown here) do not show any significant peak at the times of the jumps.Inconsistencies in the applied correction models of the ionosphere, troposphere and tides can be ruled out as the origin of these jumps.None of the time series of these features show variations exceeding a few centimetres.Future studies including additional data sets will hopefully show whether these jumps are related to remaining processing issues or physical processes.

ICESat
Prior to the validation of ICESat elevation data, we first determine the ICESat laser campaign biases as described in Sect.3.2.2.The resulting biases are given in Table 3.The simultaneously derived surface elevation change rate ḣ from Eq. ( 4) amounts to 0.0 ± 0.2 cm yr −1 .This is new, independent evidence for the stability of the surface elevation above Lake Vostok.It confirms our results in Sect.2.5 and those of previous studies (Richter et al., 2008(Richter et al., , 2014a)).It also justifies the assumption of a stable surface made by Ewert et al.
(2012) as a precondition for the campaign bias determination.It is, therefore, not surprising that our biases are very similar to the updated set of Ewert et al. ( 2012) for R33 including the Gaussian-Centroid (G-C) correction, presented by Richter et al. (2014a).The major difference is an offset of about 5 cm.This arises from the fact that in this study, we perform an absolute calibration.
The chronological sequence of the laser campaign biases implies a trend ḣb which distorts any determination of surface elevation rates if the biases are not applied.This trend over the entire ICESat operational period amounts to 1.17 ± 0.34 cm yr −1 .Table 4 and Fig. 7a compare the results of our new set of laser campaign biases with those of recent publications.For consistency, we limit this comparison to publications using either R33 data including G-C correction or R34.Hofton et al. ( 2013) used an internal crossover adjustment over the low-precipitation zone of the East Antarctic ice sheet (EAIS) as well as absolute calibration using an Figure 5. Histograms, means and standard deviations of crossover differences between different CryoSat-2 data sets and kinematic GNSS profiles for four zones of characteristic surface slope.The histogram ranges are the same as in Fig. 4 (Envisat) for comparability.The displayed crossover differences contain uncertainties in the kinematic GNSS profiles (4-9 cm) and possible elevation changes between the observation epochs of both techniques in addition to the uncertainty in the CryoSat-2 data.
ICEBridge lidar profile along the maximum latitudinal extent of the ICESat mission (86 • S).To account for any possible elevation changes, the authors apply corrections for glacial isostatic adjustment (GIA) and firn densification.It is interesting to note that, in principle, both sets of their biases are very similar to our results but their artificial trends ḣb are 0.2-0.6 cm yr −1 larger than ours (Fig. 7a, column ḣb in Table 4).This might be due to some unmodelled effects in the firn densification model and/or GIA correction of Hofton et al. (2013).Nevertheless, both of their sets of laser campaign biases (obtained in different regions of Antarctica and with different methods) agree with our results within their stated accuracies.
In contrast, the trends obtained by Urban (in Scambos and Shuman, 2016) and Zwally et al. (2015) from calibrations over the ocean differ significantly from the results of calibrations over Antarctica.Zwally et al. (2015) calculated offsets from open water and thin ice in leads and polynyas in polar sea ice and used them to determine elevation changes over Antarctica.They obtain an elevation change rate ḣ for Lake Vostok of 2.02 cm yr −1 .This contradicts the results of this study, but also those of two independent data sets in Richter et al. (2014a), i.e. static GNSS observations and kinematic GNSS profiles using snow mobiles (compare Richter et al., 2016).It is interesting to note, however, that the trends implied by our laser campaign biases and those of Zwally et al. (2015) differ by 2.11 cm yr −1 .This explains the discrepancies of the elevation change rates obtained by Zwally et al. (2015) over Lake Vostok as a result of the applied set of laser campaign biases.The choice of these biases influences the derivation of elevation change rates from ICESat over the entire Antarctic ice sheet.Hence this also explains the disparity between their ICESat-derived mass budget and the massbalance estimates of many other studies (e.g.Shepherd et al., 2012;McMillan et al., 2014;Martin-Espanol et al., 2016), especially in East Antarctica as documented in Scambos and Shuman (2016) and Richter et al. (2016).Urban et al. (2013) (as updated in Scambos and Shuman, 2016) obtained significantly smaller biases using the global ocean as a reference surface.The 2016 paper also discusses that due to different GLAS sensor saturation levels, those ocean-derived biases are not applicable over a high-albedo ice sheet surface.
After applying the laser campaign biases as corrections to the ICESat surface elevations, we calculate the crossover differences with respect to the kinematic GNSS profiles.The results in Fig. 7b show the very high accuracy of the ICESat data, even in the coastal zone.The crossover differences in the less sloping regions indicate that both data sets have practically the same precision (compare Table 1).Here, the biascorrections have only a minor influence on the standard deviations (e.g.11.1 cm for the uncorrected, 10.7 cm for the corrected data over Lake Vostok) but become much more important when the temporal distribution of the data is analysed.Close to the coast we observe a small increase in standard deviation (approx.30 cm for slopes exceeding 0.5 • ).This might be an effect of the increased surface roughness which affects the interpolation of the elevation to the crossover point.A Table 3. Intra-campaign precision RMS X from crossovers within each ICESat laser campaign and obtained campaign biases b for release 34, derived from a combined crossover adjustment of ICESat elevations and kinematic GNSS profiles in the region of subglacial Lake Vostok (to be subtracted from the elevations for correction).The true surface elevation change ḣ, estimated simultaneously in Eq. ( 4), is 0.0 ± 0.2 cm.part of this increased noise could also be explained by topographic effects within the 65 m laser footprint.

Laser
To avoid the influence of surface elevation changes, previous studies as Siegfried et al. (2011) or Kohler et al. (2013) validated only laser campaigns in very close temporal proximity (within some months) to their profiles.With a standard deviation of 10-20 cm in the majority of the area, we can confirm the results of Kohler et al. (2013).However, due to the different ICESat release versions used, their mean offwww.the-cryosphere.net/11/1111/2017/The Cryosphere, 11, 1111-1130, 2017 sets are not comparable to our results.Siegfried et al. (2011) use measurements on snowmobiles near Greenland summit for a similar validation.Their spread of ICESat -GPS differences is significantly lower than in our study.However, their assessment is limited to a single ICESat repeat pass section of 6 km length.In contrast to those studies, our profiles cover a much longer temporal range.This allows us to pinpoint the maximum range of elevation changes ḣ in this area to very tight limits.Hence, we do not correct for elevation changes in this comparison.Especially in the coastal zone, this could explain the increased offset as well as the larger noise of the crossovers.The obtained offset and standard deviation, in turn, constrain the magnitude of possible elevation changes in these zones in addition to the results in Sect.2.5.A deeper view into the temporal variability of the crossover differences in the Lake Vostok region is given in Fig. 7c.To show the variability during the different GNSS seasons, we selected a single ICESat campaign (L3C, which has the highest precision RMS X ; see Table 3) as reference and analysed the spread of the crossovers over time.The monthly averages vary by less than 10 cm and, within their standard deviations, nearly all of them can be considered to be zero.The spread of the values is very likely a result of the remaining uncertainties of the antenna height reduction in the GNSS profiles.Interannual variations in accumulation (Ekaykin et al., 2004, σ < 5 mm yr −1 ) or a water discharge from Lake Vostok (Richter et al., 2014b) can be ruled out as possible causes.Due to the long temporal base, those data sets nevertheless allow a precise trend ḣ of −0.4 ± 0.4 cm yr −1 to be derived.Together with the result of Sect.2.5 and the estimate for the surface elevation change rate ḣ, obtained from Eq. ( 4) in the laser campaign bias estimation, those are three very consistent results from different methods of obtaining ḣ from our GNSS profiles.

Data
Our kinematic GNSS profiles allow us to validate not only altimetric surface elevations but also derived products such as gridded digital elevation models (DEMs).These products are used in a wide range of applications in polar sciences.In some cases, such as the topographic correction in repeat track analysis (Moholdt et al., 2010) or the estimation of drainage basins (Zwally et al., 2012), only the precision of elevation differences between neighbouring cells is important.Other applications, such as the derivation of ice thickness at the grounding line (Rignot et al., 2008), depend on the absolute accuracy.
We validate four DEMs of Antarctica derived from satellite altimetry.The 500 m resolution DEM from data of the ICESat mission by DiMarzio et al. (2007) (further called ICESat-DEM) was a milestone for many applications.Compared to previous, SRA-based DEMs, it provided a "greater latitudinal extent and fewer slope-related effects".Nevertheless, a weak spot was the coarse cross-track spacing, especially for applications in coastal regions.The DEM provided by Bamber et al. (2009) (Bamber-DEM) overcame this problem by combining the high accuracy of ICESat with the high spatial resolution of ERS-1.The DEM produced by the Bedmap2 project (Bedmap2-DEM) combined the Bamber-DEM with regional models in the margins, the ice shelves and the Antarctic Peninsula (Fretwell et al., 2013).To make the Bedmap2-DEM comparable, we converted the elevations from the GL04C geoid to the WGS-84 ellipsoid reference surface.Even though it should be identical to the Bamber-DEM for the major part of the region, we included this model to show the loss of accuracy by rounding the elevations towards integer metres as has been done for Bedmap2.
In addition to these ICESat-dominated DEMs, Helm et al. (2014) compiled the first 3 years of the CryoSat-2 mission to a new DEM (CryoSat-DEM).With its improved design, the radar altimeter of CryoSat-2 is capable to provide very precise elevations in the margins.Furthermore, the high data density due to the orbit configuration allows for a very homogeneous data set.There is almost no need to fill data gaps due to the very small across-track spacing.

Methods
By interpolating the DEM grid to the positions of the individual GNSS measurements using bicubic interpolation, we obtain an elevation difference for every single GNSS data location.Hence, the DEM validation relies on many more elevation differences than the validation of the altimetry profiles themselves, where the heights could only be compared at crossover locations.
To facilitate the comparison with the results from Sect. 3, we subdivide these elevation differences into the same zones described in Sect.3.2.1.In the validation of DEMs, special attention has to be paid to interpolation errors.The high resolution of 500 m of the ICESat-DEM seems reasonable when working with cells which contain measurements.However, no data exist within the almost 20 km gaps between the altimeter tracks.Thus, for a closer look into these interpolation errors, we examine the dependence of the elevation differences from the distance to the nearest track.

Results
The results of the validation (Fig. 8) of the ICESat-DEM show that over Lake Vostok the accuracy is close to that of the original ICESat elevations.This indicates that in this case of exceptionally smooth topography the interpolation error is negligible.However, with increasing slope the standard deviation of this DEM grows fast due to the scarce across-track spacing of the altimetric profiles.In the Bamber-DEM, the inclusion of the additional ERS-1 data reduces this interpolation error between the ICESat tracks by 50 % in terms of standard deviation.In the coastal zone, the gain in accuracy is minor as the precision of the radar altimetry data deteriorates.
The comparison of the Bedmap2-DEM with the Bamber-DEM reveals remarkable differences.Firstly, the rounding to integer elevations in Bedmap2 increases the standard deviation by several decimetres.Secondly, a constant offset of −1 m becomes evident in the Bedmap2 data set.This seems to be an issue in the compilation procedure as the original Bamber-DEM shows a good agreement in terms of mean difference with our GNSS profiles in all regions.Finally, in the coastal zone the two models behave differently.In Bedmap2 regional elevation models have been included here, but the sparse sampling of these areas by our profiles (see Fig. 3) does not allow for a representative evaluation.
The comparison of the CryoSat-DEM with the ICESatbased models proves that SRA with the advanced instrument  design is able to provide accurate elevation information even in steep topography too.As for the CryoSat-2 altimetry data, a slope-dependent offset is observed which presumably results from the systematic underrepresentation of local depressions.
In order to shed light on the relation of the DEM accuracy with the surface slope, the median absolute deviations (DEM vs. kinematic GNSS) are binned in 0.1 • slope intervals (Fig. 9).The medians of the entire set of deviations (thick black lines) reveal a significant increase with the slope.Splitting the deviations into subsets according to the distance to the nearest ICESat track the impact of the interpolation error on a particular DEM becomes evident (panels a, b).It demonstrates the much greater dependence of the ICESat-DEM on the across-track distance compared to the Bamber-DEM.While the ICESat-DEM yields superior accuracy close to the tracks (i.e.distances < 1 km) due to its small grid interval (500 m), its median deviations exceed 10 m at large slopes (> 0.5 • ) and distances (> 6 km).
The dense spatial sampling of the CryoSat-2 altimetry data usually yields observed elevations for each cell of the CryoSat-DEM.According to the histograms obtained for the four zones, the slightly larger deviations of this DEM compared to the other models are mainly due to the slopedependent offset.Helm et al. (2014) provide an error map based on a validation with ICESat and airborne elevation data.These uncertainties are shown in Fig. 9c as a function of slope.Their good agreement with our results confirms them as realistic.

Conclusions
Between 2001 and 2015 we logged kinematic GNSS data along nine scientific convoys and derived more than 30,000 km of surface elevation profiles.We resolved the challenges of the GNSS processing, such as the very long baselines and peculiar vehicle dynamics in soft snow.Our elevation profiles have accuracies between 4 and 9 cm for a single data point.Over Lake Vostok crossover differences yield a mean elevation change rate of −0.1 ± 0.5 cm yr −1 .This confirms the results of Richter et al. (2014a) and qualifies this area as a calibration site for satellite altimetry.
A crossover analysis with three different Envisat elevation data sets reveals the impact of different processing strategies in satellite radar altimetry.Concerning the slope correction, the relocation method is clearly superior to the direct method, reducing elevation errors by about 66 % in terms of standard deviation.Threshold retrackers (ICE1/OCOG) outperform functional fit retrackers by up to 50 % in standard deviation.A similar analysis with CryoSat-2 LRM mode data confirms this finding.Hence, we recommend threshold retrackers for the determination of elevation time series because of its significant suppression of the snowpack-related pseudo elevation variations.
ICESat elevation data and our kinematic GNSS profiles are comparable in their accuracy, even close to the coast.This comparison also constrains the magnitude of temporal elevation changes.A combined crossover adjustment above Lake Vostok yields a new set of ICESat laser campaign biases that no longer depends on an a priori assumption of a stable surface elevation.The surface elevation change rate of 0.0 ± 0.2 cm yr −1 , obtained here, nevertheless proves that this assumption is correct to a very high level of certainty.The correction of ICESat elevation data for the laser campaign biases is crucial for estimates of surface elevation change rates and the according mass balance.The discrepancy of the bias-induced trend between Zwally et al. (2015) and our result as well as other recent studies (Hofton et al., 2013;Richter et al., 2014a) is approximately −2 cm yr −1 .This is very likely the cause for the significantly more positive Antarctic ice mass balance estimate from ICESat obtained by Zwally et al. (2015) compared to other authors (see Scambos and Shuman, 2016).
The validation of four digital elevation models demonstrates the reduction of interpolation errors achieved by Bamber et al. (2009) by complementing ICESat elevations with radar altimetry data.The advanced instrument design and high spatial resolution of CryoSat-2 now permits radar altimetry to provide DEMs (Helm et al., 2014) of similar accuracy, avoiding extensive interpolation.
Data availability.The full set of surface elevation profiles obtained from kinematic GNSS measurements is available for download on the data server PANGAEA (https://doi.org/10.1594/PANGAEA.869761).
Competing interests.The authors declare that they have no conflict of interest.
Figure 1.(a) Overview of the kinematic GNSS profiles (colour coded by their chronological sequence) and the GNSS reference stations used in the differential processing.(b) Detailed map of the profiles in the area of subglacial Lake Vostok (outline in grey, hydrostatic equilibrium area in black).(c) Convoy vehicle of type STT-2 Kharkovchanka-2 (profile K10B) with antenna mounted on the container above the cabin (red circle).(d) Convoy vehicles Kässbohrer PistenBully (profiles K14A and K14B) with antenna mounted on top of the cabins.

Figure 2
Figure 2. (a) Elevation change rate ḣ from crossovers between different seasons (20 km block mean values).(b) Maximum observation time span of crossover differences within each block.(c) ḣ in the Lake Vostok region.Crossovers on the convoy track or with t < 5 yr are considered uncertain and thus plotted half sized.(d) Distribution of the valid crossovers ( t>=5 yr, no anthropogenic surface disturbances) in the Lake Vostok region (100-108.5• E, 76-79 • S).

Figure 3 .
Figure3.Crossover differences between Envisat (ESA, relocated, ICE1 retracker) and the kinematic GNSS profiles.The black contour lines mark the borders between the different zones of slope (0.15 and 0.5 • ) and the hydrostatic equilibrium area of Lake Vostok (outline in grey fromPopov and Chernoglazov, 2011).

Figure 6 .
Figure 6.(a-d) Monthly averages of crossover differences between different versions of CryoSat-2 data and kinematic GNSS profiles (black) and the corresponding backscatter σ 0 (red) within the hydrostatic equilibrium area of Lake Vostok.(a) Baseline B data set, (b-d) Baseline C using the three different retrackers applied.The box at the bottom of each plot gives the overall elevation trend.

Figure 7 .
Figure 7. (a) ICESat laser campaign biases (release 34) determined in the present study in the Lake Vostok region (100-108.5• E, 76-79 • S) from inter-campaign crossovers and from crossovers with kinematic GNSS profiles (in dark blue) compared to other recently published sets of biases and their trends.Error bars represent the standard deviation (if given) for the respective data set.All bias sets have been reduced by their mean value, given in the legend.(b) Histograms, means and standard deviations of crossover differences between ICESat and the kinematic GNSS profiles for four zones of characteristic surface slope.The range of the histograms is kept fixed to ±1 m.The displayed crossover differences contain uncertainties in the kinematic GNSS profiles (4-9 cm) and possible elevation changes between the observation epochs of both techniques in addition to the uncertainty in the ICESat data.(c) Monthly averages of crossover differences between ICESat campaign L3C and the kinematic GNSS profiles in the region used for the bias determination.

Figure 8 .
Figure 8. Histograms, means and standard deviations of differences between four different DEMs and kinematic GNSS profiles for four zones of characteristic surface slope.

Figure 9 .
Figure 9. Median of absolute differences between various DEMs and the kinematic GNSS profiles binned by surface slope (thick black line).(a-b) Additionally for DEMs based on ICESatcoloured lines show the differences for several ranges of distances to the laser altimeter tracks (colour coded).(c) For the CryoSat-DEM the magenta line shows the estimated error given by the uncertainty map which comes with the DEM.
Table 1 contains detailed information about each individual profile.

Table 1 .
Overview of the kinematic GNSS profiles.Estimates of precision consist of the mean elevation error from GNSS processing (from differences between combined and single baseline solutions) RMS BL and the obtained mean formal precision of snow surface elevation RMS S .The empirical mean crossover difference between independent profiles of one season RMS X gives an estimate for the accuracy.Averages and (in case of more than one) the standard deviation of the crossover differences with Operation ICEBridge ATM at 26 November 2013 is given in the last column.Duration dates are given as yyyy-mm-dd. correspond

Table 4 .
Zwally et al. (2015)nferred from different sets of ICESat laser campaign biases weighted according to their standard deviation (if given).Results for Hofton (2013) differ from the originally given values as these authors applied unit weights.The trends in the second column have been calculated using only the laser campaignsZwally et al. (2015)employed for their study.