Improved retrieval of land ice topography from CryoSat-2 data and its impact for volume-change estimation of the Greenland Ice Sheet

. A new methodology for retrieval of glacier and ice sheet elevations and elevation changes from CryoSat-2 data is presented. Surface elevations and elevation changes determined using this approach show signiﬁcant improvements over ESA’s publicly available CryoSat-2 elevation product (L2 Baseline-B). The results are compared to near-coincident airborne laser altimetry from NASA’s Operation IceBridge and seasonal height amplitudes from the Ice, Cloud, and Elevation Satellite (ICESat). Applying this methodology to CryoSat-2 data collected in interferometric synthetic aperture mode (SIN) over the high-relief regions of the Greenland Ice Sheet we ﬁnd an improvement in the root-mean-square error (RMSE) of 27 and 40 % compared to ESA’s L2 product in the derived elevation and elevation changes, respectively. In the interior part of the ice sheet, where CryoSat-2 operates in low-resolution mode (LRM), we ﬁnd an improvement in the RMSE of 68 and 55 % in the derived elevation and elevation changes, re-spectively. There is also an 86 % improvement in the magnitude of the seasonal amplitudes when compared to amplitudes derived from ICESat data. These results to in dielectric

Abstract. A new methodology for retrieval of glacier and ice sheet elevations and elevation changes from CryoSat-2 data is presented. Surface elevations and elevation changes determined using this approach show significant improvements over ESA's publicly available CryoSat-2 elevation product (L2 Baseline-B). The results are compared to near-coincident airborne laser altimetry from NASA's Operation IceBridge and seasonal height amplitudes from the Ice, Cloud, and Elevation Satellite (ICESat).
Applying this methodology to CryoSat-2 data collected in interferometric synthetic aperture mode (SIN) over the highrelief regions of the Greenland Ice Sheet we find an improvement in the root-mean-square error (RMSE) of 27 and 40 % compared to ESA's L2 product in the derived elevation and elevation changes, respectively. In the interior part of the ice sheet, where CryoSat-2 operates in low-resolution mode (LRM), we find an improvement in the RMSE of 68 and 55 % in the derived elevation and elevation changes, respectively. There is also an 86 % improvement in the magnitude of the seasonal amplitudes when compared to amplitudes derived from ICESat data. These results indicate that the new methodology provides improved tracking of the snow/ice surface with lower sensitivity to changes in nearsurface dielectric properties.
To demonstrate the utility of the new processing methodology we produce elevations, elevation changes, and total volume changes from CryoSat-2 data for the Greenland Ice Sheet during the period January 2011 to January 2015. We find that the Greenland Ice Sheet decreased in volume at a rate of 289 ± 20 km 3 a −1 , with high interannual variability and spatial heterogeneity in rates of loss. This rate is

Introduction
In April 2010, the European Space Agency (ESA) launched CryoSat-2 tasked with monitoring the changes of the Earth's land and sea ice. CryoSat-2 carries a new type of delay/Doppler radar altimeter (Raney, 1998) referred to as SIRAL (SAR Interferometric Radar Altimeter). SIRAL operates in two different modes over land ice. Over the interior part of the ice sheets it operates as a conventional pulse limited radar system, referred to as the low-resolution mode (LRM). In more complex high-sloping terrain the system uses a second antenna to operate in interferometric synthetic aperture radar (SIN) mode. This new feature allows the satellite to monitor changes in complex terrain, including ice caps, glaciers, and the high-relief marginal areas of the ice sheets. Such areas are sensitive to changes in climate and contribute greatly to current rates of sea level rise (e.g., Gardner et al., 2013;Shepherd et al., 2012).
Ku-band radar altimeters are insensitive to cloud cover providing superior coverage compared to laser altimeters (e.g., ICESat) but experience significant amounts of volume scattering, which is controlled by the time-evolving dielectric properties of the near-surface snow, firn, and ice (Lacroix et al., 2008;Remy et al., 2012). These effects can have large implications for the determination of mass change over a wide range of both spatial and temporal scales. Changing snow Published by Copernicus Publications on behalf of the European Geosciences Union.
conditions can introduce time-varying biases in the data that, in combination with the radar signals interaction with the surface, introduce large elevation biases (0.5-1 m) (Nilsson et al., 2015a). This, combined with other factors such as processing methodology and surface topography, makes it difficult to measure small changes for much of the word's icecovered regions (Arthern et al., 2001;Gray et al., 2015;Nilsson et al., 2015b).
The mitigation of these effects in the processing of radar altimetry data is required for improved accuracy of derived temporal and spatial changes in surface elevation of glaciers and ice sheets. Several studies have proposed different approaches to assess these effects and improve the retrieval process of surface elevation and elevation changes from radar altimetry data. These include different approaches to waveform retracking (Davis, 1993(Davis, , 1997Gray et al., 2015;Helm et al., 2014) and empirical corrections to the estimated surface elevation changes (Davis and Ferguson, 2004;Flament and Rémy, 2012;Sørensen et al., 2015;Wingham et al., 2006b;Zwally et al., 2005Zwally et al., , 2011. Relatively little work has been done to assess methods for improving elevation and elevation changes derived from ESA's CryoSat-2 data (Abulaitijiang et al., 2015;Gray et al., 2013Gray et al., , 2015Helm et al., 2014).
Here we conduct a thorough analysis of CryoSat-2 SIN and LRM waveform retracking, filtering and processing methodologies to design an optimal processing methodology for CryoSat-2 elevation retrieval over both smooth and complex ice-covered terrain. We then analyze two different approaches in determining surface elevation and volume changes from the scattered CryoSat-2 elevation retrievals. The overarching goal of this work is to develop robust and accurate elevation retrieval algorithms that are less sensitive to changes in surface and subsurface scattering properties.
The new processing scheme is applied to estimate elevation and volume changes of the Greenland Ice Sheet for the period January 2011 to January 2015, using two independent methods to characterize the differences of the results due to methodology. The results are compared to change estimates obtained from the ESA L2 Baseline-B surface elevation product (Bouzinac et al., 2014), high-accuracy airborne data from NASA IceBridge Airborne Topographic Mapper, and seasonal height amplitudes estimated from Ice, Cloud, and Elevation Satellite (ICESat) data.

Low-resolution mode (LRM)
The LRM data is used over the interior parts of the ice sheet, which mostly consist of low-sloping terrain. Here, SIRAL operates as a conventional pulse-limited radar system with a transmission frequency of 13.6 GHz (Ku-band) and has a pulse-limited footprint (PLF) radius of approximately 1.5 km and a beam-limited footprint (BLF) radius of approximately 7.5 km over flat terrain (Bouzinac, 2014). The gentle terrain allows for accurate mapping of the surface elevation of the ice sheet down to decimeter-level (Brenner et al., 2007). Within the LRM waveform we define the location of the surface from the leading edge of the waveform based on a fraction of the maximum amplitude of the received power. This approach is commonly referred to as a threshold retracker. Following Davis (1997) we use a 20 % threshold to define the location of the surface. Davis (1997) argued that a 20 % threshold represents the best compromise between waveforms that are entirely dominated by either volume or surface scattering, making it suitable for obtaining estimates of surface elevation for most parts of the Greenland Ice Sheet.
The CryoSat-2 LRM radar waveforms suffer from measurement noise, in the form of speckle noise. Furthermore, over the steeper parts of the LRM area the range gate tracking loop can lose track of the surface, producing unusable waveforms. To remove bad or loss-of-track waveforms the radar waveform (20 Hz) is first filtered using a zero-phase low-pass filter to reduce speckle noise on a line-by-line basis. The signal-to-noise ratio (SNR) of the waveform is then estimated and if the SNR < 0.5 dB the waveform is rejected. The SNR threshold was empirically chosen to obtain a good trade-off between the quality of the measurements and sampling.
Before the waveform can be retracked the first surface return (first major peak) is identified within the range gate window. A copy of the waveform is heavily smoothed to remove small-scale surface roughness signals, keeping the overall surface signal intact. The range gate index of the first peak from the copy is then used to extract the leading edge of the original low-pass filtered waveform. Only leading edges with a peak index above 20 are used in the retracking, as peaks before that can indicate troublesome surface ranging. The extracted leading edge is then oversampled by a factor of 100 (cf. Gray et al., 2013;Helm et al., 2014), and the range R between the surface and satellite is determined based on the 20 % threshold computed according to Davis (1997). The range is then corrected for several atmospheric and geophysical effects relevant to land ice studies according to Bouzinac (2014). The surface elevation H of the topography, relative to the WGS84 ellipsoid, is estimated as H = A − R, where A is the altitude of the satellite.
The measured surface return over a sloping surface does not originate from the satellites' nadir location, but from the "point of closest approach" (POCA) to the spacecraft (Brenner et al., 1983). These off-nadir returns can introduce a large-range bias to the surface, depending on the magnitude of surface slope, ranging from 0 to 120 m (Brenner et al., 1983) as the measured surface height is mapped to an erroneous position (i.e., the nadir position). To mitigate the effect of this error we correct the measured range and location to the POCA point using an a priori DEM, following the approach of Bamber (1994). In contrast to the previous study we account also for the local surface curvature, as The Cryosphere, 10, [2953][2954][2955][2956][2957][2958][2959][2960][2961][2962][2963][2964][2965][2966][2967][2968][2969]2016 www.the-cryosphere.net/10/2953/2016/ Remy et al. (1989) showed that accounting for surface curvature in addition to surface slope significantly improves results. The surface slope, aspect, and curvature are estimated from an a priori DEM. The GIMP elevation model (Howat et al., 2014) was used to derive surface parameters for the slope-induced error correction in the LRM data. The DEM was resampled to 2 km resolution, using bilinear interpolation prior to parameter estimation, which provided the lowest root-mean-square error and further corresponds to the pulselimited footprint of the LRM data.

Interferometric synthetic aperture radar (SIN) mode
The SIN mode is used over the marginal areas of the ice sheets and other smaller glacier-iced areas. In these areas the SIRAL altimeter operates as a delay/Doppler radar system (Raney, 1998). The delay/Doppler radar allows for higher along-track resolution compared to conventional altimetry, resulting in 350 m resolution in along track and 1500 m across track. In ordinary SAR operation only the amplitude of the radar echo is measured and the phase content is discarded or ignored. With the inclusion of a second antenna on CryoSat-2, interferometric SAR can be performed. The difference in path length between the POCA and the individual antennas introduces a phase shift between the two retrieved signals that can be related to the angle of arrival (look angle). The look angle can in turn be used to resolve the across-track (across-antenna) location of the echo. Multi-look processing is applied to ESA's L1B waveform product (Bouzinac, 2014) to reduce the noise in the SIN waveform but it is still affected by speckle noise, as is the case for the LRM waveforms. To mitigate this effect, and to help identify the leading edge of the first return, we apply speckle reduction filtering and leading edge extraction of the SIN waveforms in the same way as for the LRM data processing, with minor changes due to differences in range gate resolution. In this case, compared to the LRM retracking algorithm, only leading edges with a peak index in the range of 100-350 are used for retracking the radar waveform (P w ).
The estimated coherence C of the multi-looked waveforms is then filtered in two stages: (i) all coherence measurements larger than one are set to zeros (coherence values larger than one exists in the L1B product due to unknown reasons); (ii) the coherence array, as a function of range, is filtered using a 2-D 5×5 Wiener filter to remove high-frequency noise. The filtering of the waveform and the coherence is applied to remove noise in the recreation of the interferogram, discussed below.
The measured differential phase φ of the return signal is affected by phase ambiguities, for example, a sudden shift of 2π in the measured phase. To reduce phase noise and aid the phase, an unwrapping of the radar interferogram I is per-formed according to Gray et al. (2013): (1) The interferogram is then filtered using a wavelet-based denoising technique, where the real and imaginary parts of the interferogram are filtered separately. The unwrapping of the interferogram allows for indirect filtering of the phase without being affected by the phase ambiguities. Phase filtering is an important consideration as it has a direct effect on accuracy of the position of the ground echo. We selected a bi-orthogonal wavelet as the mother wavelet to produce the wavelet coefficients decomposed into three levels.
Soft thresholding was applied to detail coefficients, using a heuristic threshold rule to remove noise at every level. This was done on a line-by-line basis. The final filtered differential phase was then recovered by where I Re f is the real part of the interferogram and I Im f is the imaginary part. To resolve the phase ambiguities the filtered phase measurements require unwrapping. The phase unwrapping is done on a line-by-line basis in two directions starting from the center of gravity of the waveform (Wingham et al., 1986).
The return power distribution of a delay/Doppler radar system shows an important distinction from those from conventional pulse-limited radar systems. Here, the point corresponding to the mean surface is not located at the halfpower point on the leading edge, but rather closer to the maximum (Wingham et al., 2006a). Therefore, a new retracker has been developed, closely related to the one used in Gray et al. (2013), to allow for adaptive retracking of the upper parts of the leading edge of the SAR waveform. The algorithm follows the main concept of the threshold retracker, developed by Davis (1997), but instead of a pre-defined threshold it tracks the maximum gradient of the leading edge of the oversampled waveform. We refer to this approach as the "leading-edge maximum gradient retracker" (LMG).
The surface returns are geolocated using the across-track look-angle θ estimated from the differential phase at the retracking point according to Wingham et al. (2006a). This, in combination with the viewing geometry, is used to define the location of the surface return on the ground using basic across-track interferometric principles. We correct θ for the interferometer surface slope error by applying the look-angle scaling factor estimated in Galin et al. (2013) and for the platform roll angle.
The along-track differential phase estimate, interpolated to the retracking point, is affected by phase ambiguities not corrected for during the phase unwrapping procedure. To reduce residual phase ambiguities an a priori DEM (GIMP) is used to extract the DEM surface elevations at the nadir, resampled to 500 m resolution (corresponding roughly to the alongtrack sampling), and echolocation using bilinear interpolation. Over a sloping surface the surface return should always www.the-cryosphere.net/10/2953/2016/ The Cryosphere, 10, 2953-2969, 2016 come from a position upslope from the nadir point. Therefore the following relation must hold where (H echo > H nadir ) or for a more practical application (H echo − H nadir ) > ε DEM , where ε DEM is the uncertainty of the DEM used. If this relation is violated 2π is added or subtracted to the individual along-track phase estimate, depending on the sign. A final step is applied to correct for any lingering phase ambiguities not corrected by the a priori DEM. This step uses the assumption that the along-track phase should follow a consistent pattern over most of the satellite ground track. Hence, any large discrepancies from the overall pattern of the along-track phase would indicate an ambiguity. The ambiguity is detected by computing the residuals of the alongtrack phase by removing a smoothed version of the differential phase. If any of the residuals have a magnitude larger than π they are considered ambiguous and thus corrected by adding or subtracting 2π .

Surface-fit method
The surface-fitting method (SF) is based on fitting a linear model to the elevations as a function of time and space inside a search radius of 1 km (e.g., Howat et al., 2008;Moholdt et al., 2010;Sørensen et al., 2011;Wouters et al., 2015). The linear model consists of a time-invariant (static) bi-quadratic surface model to account for variable topography inside the search radius and a time-variant part used to extract the temporal change in elevation. The model consists of a total of seven parameters, where six of the parameters (a coefficients) describe the bi-quadratic surface-modeling function, ∂h/∂t the linear elevation change rate, t time in decimal years, t 0 the mean time inside the footprint, and ε the residuals from the linear regression.
h (x, y, t) = a 0 + a 1 x + a 2 y + a 3 xy + a 4 x 2 + a 5 y 2 The algorithm estimates the elevation change at every echolocation (or grid node if desired) in the data set. In each solution the signal amplitude and phase are also estimated by fitting a seasonal signal model to the surface-fit elevation residuals, according to where h is the elevation residuals estimated from the surface-fit model, s 0, 1 are the model coefficients, and t is the time. The amplitude A is then defined as A = s 2 0 + s 2 1 and the phase P as P = tan −1 s 1 s 0 . To remove outliers an iterative 3σ filter is used in the full model solution i.e., the topography, trend, and seasonal signal are removed, using a maximum of five iterations. For each iteration residuals (full-model) with an absolute value larger than 10 m are removed, as seasonal changes larger than 10 m are not expected (Moholdt et al., 2010;Qi and Braun, 2013). The data inside the 1 km cap are weighted according to their distance from the estimation point according to where W is the estimated weight, d the distance, and ρ the correlation or resolution parameter set to 500 m. The weighting allows the solution to better reflect local signal dynamics at the prediction point. Local elevation time series are further computed from the elevation residuals and elevation trend for each solution according to where t is the time epochs inside the search cap, t 0 is the mean time of t, ∂h/∂t is the estimated elevation change rate, and ε is the elevation residual at each time epoch. The elevation changes estimated from the surface-fitting method are then culled to remove outliers before spatial gridding. Elevation changes with an absolute trend larger than 15 m a −1 are removed. The resulting surface elevations are binned at 5 km resolution for outlier-editing purposes. For each cell the local spatial trend is modeled as a bilinear surface, and removed. The residuals are then edited using an iterative 3σ filter until the RMS (root mean square error) converges to 2 %.

Crossover method
The crossover method is used to derive the surface elevations at the intersection point between an ascending and descending satellite ground track separated in time (Brenner et al., 2007;Khvorostovsky, 2012;Zwally et al., 1989). The surface elevations and times are then estimated at the crossover location for each track by linear interpolation of the two closest data points for each ascending and descending track. The crossover height difference is then estimated by taking the height difference between the two tracks according to where h 1 and h 2 are the surface heights at the crossover location at time epoch t 1 and t 2 , respectively, and ε is the random measurement error, including orbital, range, and retracking errors. This approach produces crossover height differences with scattered time epochs ranging from 0 to 4 years. CryoSat-2 has a 369-day repeat orbit configuration with a 30-day subcycle, meaning that each crossover location will be revisited every 369 days and the surrounding area every 30 days.
The Cryosphere, 10, 2953-2969, 2016 www.the-cryosphere.net/10/2953/2016/ This produces annual and sub-annual crossover difference around each crossover location. This fact is used to produce elevation change rates by incorporating all multi-temporal crossover difference within a neighborhood of 2.5 km around each crossover location. The elevation change is then estimated using the same procedure described for the surface-fit method, except that a bilinear model is used to remove any spatial trends in the topography of the crossover elevations according to where dh is the crossover height difference, t the time, t 0 the mean reference time inside the footprint, a 1 and a 2 the across and along-track slope, and ∂h/∂t the elevation change rate. This produces elevation changes comparable in time and in spatial coverage with the surface-fit method. The same outlier editing schemes are applied to the crossover elevation change rates as for the surface-fit method.

Gridding of sparse elevation and elevation-change data
The gridding is done in a polar-stereographic projection with a latitude of origin at 70 • N, central longitude of 45 • W, and origin at the North Pole, and is referenced against the WGS-84 ellipsoid. The observations derived from the surface-fit and crossover methods are gridded at a resolution of 1 km × 1 km, due to the high spatial sampling. The method of least squares collocation (LSC), described in Herzfeld (1992), is used to grid the observations onto a regular grid. LSC is similar to Kriging and allows for optimal interpolation and merging of data with different accuracies, using their inherent covariance structure. The LSC algorithm uses the 25 closest data points in 8 quadrants surrounding the prediction point to reduce spatial biasing. The prediction equation consists of two terms where the first term is the actual prediction term and the second term accounts for the non-stationary part of the data, as described bŷ where C sz is the cross-covariance, C zz is the auto-covariance, N the diagonal noise matrix consisting of the a priori RMSE, and m(z) is the median value of the observations inside the search neighborhood.
The covariance of the data inside the local neighborhood is modeled as a function of distance away from the prediction point using a third-order Gauss-Markov model described below.
where r is the separation distance, C 0 is the local data variance, and α is a scaling factor estimated from the correlation length.
LSC interpolation provides a RMSE for each prediction point estimated from the modeled covariance of the data according to where the RMSE of the prediction equals σŝ = (Cŝ) 1/2 and C T sz is the transposed cross-covariance matrix. The elevation changes estimated from the surface-fit and crossover methods are interpolated using their a priori error estimated from the regression. To avoid unrealistically small errors, common in the regression errors estimated over flat terrain, a minimum error threshold is applied. Error values smaller than a specific threshold are set to the threshold value. The threshold value is representative of the overall precision of the elevation changes over flat terrain and is set to 0.2 m a −1 . The data are then gridded using a 75 km correlation length determined from the comparison of CryoSat-2 elevation to airborne measurements (Sect. 5).
The LSC algorithm is also used to generate a DEM based on the surface elevations generated from the surface-fitting algorithm. The surface elevations generated from the surface fit were used as input for the gridding algorithm. The use of surface elevations from the surface fit provides several advantages compared to the raw observations because they provide an almost equal number of observations as the raw data, have been screened for gross outliers, have been lowpass filtered using the 1 km search radius, and all reference roughly the same time epoch. Furthermore, the RMSE generated from the surface-fit-estimated surface height can be used as an a priori error for the LSC gridding procedure.
The DEM is generated using the same approach as for the surface elevation changes, as described previously in this section, with a correlation length of 25 km. Before the gridding procedure is applied elevations H < 0 and H > 3350 m are removed from the data set. Furthermore, elevations with a standard error larger than 30 m are also removed. The elevations are spatially binned into a resolution of 1000 m. Inside each cell the local surface trend is removed by fitting a planar surface, and an iterative 3σ filter is applied to the residuals to remove outliers.

Surface elevations and elevation changes from ICESat
To assess basin-scale patterns of elevation change we compare elevation changes from CryoSat-2 data to elevation changes derived from Ice, Cloud, and Elevation Satellite (ICESat) data. Here we use Release 33 (GLA06) data collected over the 2003-2009 period. The ICESat surface heights were used to generate surface elevation changes and seasonal parameters according to method M3 in Sørensen et al. (2011). The derived elevation changes were corrected for the Gaussian centroid offset (Borsa et al., 2014). Valid elevation retrievals were selected according to Nilsson et www.the-cryosphere.net/10/2953/2016/ The Cryosphere, 10, 2953-2969, 2016 al. (2015b). The ICESat elevation, seasonal amplitude, and phase are then used for comparison with CryoSat-2 and to build continuous time series using the surface-fit method described in Sect. 3.1. For the purpose of this study no correction for the inter-campaign bias was applied, as this is still an active area of investigation.

Validation
Elevation and elevation change results were generated for the entire Greenland Ice Sheet using CryoSat-2 data collected between January 2011 and January 2015 using the methodology presented in Sects. 2 and 3 (JPL product) and by applying the methods of Sect. 3 to ESA's CryoSat-2 L2 elevation products (ESA product). Surface elevations and elevation changes were validated against airborne data sets collected by NASA's Operation Ice-Bridge Airborne Topographic Mapper (ATM), obtained from the National Snow and Ice Data Center (NSIDC) in the form of the ILATM2 product. The generated elevation product has a resolution of 80 m, with a 40 m spacing along track. This mission produces both elevation and elevation changes with reported vertical accuracy of ∼ 10 cm and temporal accuracy at the centimeter level (Krabill et al., 2002). The derived surface elevations from CryoSat-2 are differenced against ATM surface elevations within 50 m of each ATM location. One month of CryoSat-2 data consistent in time with the ATM elevations is used for the validation to avoid biases due to temporal sampling and to obtain sufficient sample size. A total of 4 years of campaign data is used for the validation of the surface elevations (2011)(2012)(2013)(2014). The residuals are edited using an iterative 3σ filter to remove outliers. The accuracy and precision is estimated as the mean and standard deviation of the differences, respectively. The residual distribution is further binned according to surface slope estimated from the GIMP DEM (Howat et al., 2014) resampled to 500 m. The sensitivity to surface slope (slope error) can be identified in the standard deviation of the binned residuals and can be used to judge the quality of the produced surface elevation and elevation changes, while the binned average for the elevations can be used to determine radarsignal-penetration depth.
Surface-elevation-change rates estimated from three different time periods (2012-2014, 2011-2013, and 2011-2014) of overlapping ATM observations (Krabill, 2016) are used to validate the surface elevation changes estimated from the CryoSat-2 data. The same validation methodology applied to surface elevations is applied to surface elevation changes, with a few minor modifications. First the search radius is increased to 175 m to make it conform to the ATM elevation-change resolution of 250 m, as this search radius encloses the entire ATM grid cell. Secondly the estimated mean and standard deviation are multiplied with the individual time intervals of the validation data sets to make the errors comparable, as they differ in time span.
For the surface-fit and crossover methods, near-coincident elevation change rates were compared with ATM rates (e.g., April 2011 to April 2014). This provided three validation data sets for the surface-fit method due to its high spatial coverage. However, only the 2011-2014 validation data set could be used for the crossover method due to the lower spatial sampling of the crossovers.
The overall accuracy and precision for both the surface elevation and elevation changes are then estimated by taking the weighted mean, using the number of observations as weights, for each data set giving an average error for each For the estimated surface elevation changes from the two independent methods we find the lowest RMSE for the surface-fit method, followed by the crossover method. This differs from the findings of Moholdt et al. (2010), who found lower intrinsic errors for the crossover method compared to the surface-fit method when applied to ICESat data. The larger search radius used for our application of the crossover method most likely explains the difference in findings between the two studies. Furthermore, we find that the surfacefit method provides the largest reduction in RMSE for the JPL product, corresponding to 40 and 55 % for the SIN and LRM data, respectively.
The correlation length used to derive the number of uncorrelated grid cells, which is used to estimate the standard error, was determined from a semi-variogram analysis of the elevation-change residuals from CryoSat-2 minus ATM using the data from the surface-fit method. The comparison was done for each mode separately for all the individual campaigns and multiplied with their individual time spans. The semi-variogram was then computed from all the timeinvariant residuals for each mode to maximize the spatial coverage. Analysis of the semi-variogram showed approximate correlation lengths of 100 and 75 km for the SIN and LRM data, respectively. These correlation lengths are inside the range of those found by Sørensen et al. (2011) for their analysis of ICESat data, which was found to be between 50 and 150 km.
The main goal of this study is not to derive or compare different types of DEMs. However, to gain insight into the overall quality of our CryoSat-2-derived DEM (referred to as JPL) we compare it to three other DEMs derived from other data sets. First, we compare it to a DEM derived from ESA CryoSat-2 L2 data (referred to as ESA) gridded in the same manner as our DEM (Sect. 3.3). Second we compare it to a DEM from Helm et al. (2014), also based on CryoSat-2 data from 2011 to 2014 (referred to as AWI). Third, we compare it to a DEM from Howat et al. (2014) (which was used to derive topographical parameters and corrections for the JPL CryoSat-2 data) based on photogrammetry and altimetry data from the mid 1990s to 2010 (depending on data source) co-registered to ICESat elevation data from 2003 to 2009 (referred to as GIMP).
These data sets were then compared to IceBridge ATM elevations spanning the four different campaigns previously used for validation of the CryoSat-2 elevations. The DEM elevation was estimated at each ATM location using bilinear interpolation, and the elevation difference was computed as (DEM − ATM). No attempt was made to account for differences in DEM and ATM epochs. The estimation of the errors of the DEM was determined in the same way as for the individual CryoSat-2 surface heights. The results of the com- parison are summarized in Table 3 as the weighted average of the different campaigns. The values from each individual campaign can be found in the Supplement.
Analyzing the overall RMSE we find that the AWI produces the lowest RMSE, followed by JPL, ESA, and GIMP, due to AWI's lower standard deviation. However, the best accuracy is obtained by the JPL DEM, which shows the lowest elevation bias of all DEMs. The ESA-derived DEM shows a slightly better standard deviation than the JPL DEM, which can probably be explained by higher data density in the marginal areas for the ESA data. The difference in density is probably due to the SNR rejection criterion applied in our elevation processing. The lower standard deviation in the AWI product is most likely due to the use of lower resolution topography in many of the high-relief areas in the 1 km elevation model, producing a smoother estimate of the surface. The GIMP data set showed higher degrees of impulse noise than the other products, explaining the higher observed standard deviation. This impulse noise is attributed to the local elevation change rate, which was not accounted for in the creation of the DEM (Howat et al., 2014). Overall we find that the JPL DEM provides a suitable compromise between resolving local detail and minimizing bias. Furthermore, modification to the SNR filtering criteria will likely lead to additional improvements in the DEM.
To determine the effect of retracking on the accuracy and precision of the measured surface heights from CryoSat-2 several tests were performed over different parts of Greenland for both modes. Following the approach of Davis (1997) the accuracy (mean) and precision (standard deviation) were computed as a function of leading-edge threshold (in percent). This computation was performed using a standard leading-edge threshold retracker, referred to from now on as LTH, for both the LRM and SIN data independently. The validation was performed in the same manner as described in Sect. 5, where ATM elevations from 2013 were used as the surface reference.
For the LRM, data from April 2013 from the northern parts of Greenland, spanning the region 75-81 • N and 54-44 • W around the North Greenland Eemian Ice Drilling camp (NEEM), were used to calculate height residuals for the difwww.the-cryosphere.net/10/2953/2016/ The Cryosphere, 10, 2953-2969, 2016 ferent thresholds. This produced approximately 1000 comparison locations, which were used to calculate statistics. The same procedure was performed over Jakobshavn Isbrae, using the same time span, to calculate statistics for the SIN mode, providing roughly 2500 comparison locations. The results of this analysis, summarized in Fig. 2, show that for the LRM data the precision (as a function of threshold) follows the same behavior as observed by Davis (1997), with a decrease of precision following increasing retracking threshold. However, the most notable finding was the observed inverse relationship in precision for the SIN mode data compared to the LRM data. For the LTH algorithm in the SIN mode, we observe a clear increase in precision as the retracking threshold increases, seen in Fig. 2, stabilizing around 30-40 %.
Analyzing the accuracy derived from the different thresholds, a clear difference in apparent penetration depth of the radar signal can be observed for the two modes. For the SIN mode, below 40 % a positive bias is observed, indicating that retracker produces elevations larger than the corresponding airborne-measured heights. For thresholds larger than 40 % surface penetration of the signal is observed, which is in general closer to the surface compared to the LRM data. We attribute this to differences in the near-surface density structure covered by the two modes.
In general we conclude that for the LRM data, applying low retracking thresholds (0-30 %) reduces the magnitude of the apparent surface penetration bias and provides higher precision compared to higher thresholds. Therefore, a threshold of ∼ 20 % of the leading edge is suggested for retracking surface elevations for the LRM data, which was also previously suggested by Davis (1997) and Helm et al. (2014). However, for the SIN mode a threshold below 40 % is not recommended, as this produces a clear positive elevation bias and poor precision, as seen in Fig. 2. Analyzing the difference between the LTH and the adaptive LMG algorithm, used in the SIN mode, we find that the LMG algorithm produces superior results in precision compared to the standard LTH algorithm. Comparing the adaptive solution from LMG to the optimum threshold (40 %) found by the LTH algorithm, we find a comparable magnitude of the elevation bias and a 32 % improvement in precision with an overall 27 % reduction in RMSE using the LMG retracker. From this comparison between the two retracker algorithms we recommend the use of the adaptive threshold approach (LMG), as it produces an elevation repeatability that exceeds that of the standard threshold retracker (LTH) and provides a low penetration bias.
A case study was also performed to determine the different effects of the processing steps on the quality of the retrieved elevations. For this purpose the Barnes Ice Cap, on Baffin Island in the Canadian Arctic, was chosen due to its small size, excellent validation coverage, and because it consists mostly of super-imposed ice (reducing radar signal penetration). The ice cap saw a major IceBridge ATM campaign in 2011 providing a large number of flight tracks (spanning in both north-south and east-west directions) suitable for validating CryoSat-2 data. The result of this case study, which is detailed in the supplementary material (i.e., Table S1 in Supplement) shows that the filtering of the differential phase has the highest impact on the overall accuracy of the observation, reducing the RMSE by 12 %, followed by the ambiguity correction. This shows the importance of these steps, as they can have important implications for the overall quality of the retrieved elevations. This is especially true in high-relief areas where small changes in the look angle, or an introduced phase ambiguity, can produce large elevation errors ranging from 0 to 100 m in elevation (Brenner et al., 1983).

Error analysis
To compute volume change errors for the two methods we divide the error budget into two main components: (1) the observational standard error (ε obs ) and (2) the interpolation standard error (ε int ).
The observational error budget is estimated using the root-mean-square error (RMSE) of the difference between CryoSat-2 and airborne elevation change differences, as described in Sect. 5. The RMSE is estimated separately from the two different modes, with the total volume change error (ε vol ) being computed as the root sum of squares (RSS) of the standard elevation change error of the two modes and their corresponding area, according to where A lrm and A sin are the areas covered by each mode. The ε lrm and ε sin are the standard elevation change errors of the LRM and SIN data computed from the airborne validation data sets. The observational elevation-change error is estimated from the residual elevation change differences in Table 2 for the two methods. The RMSE from the LRM/SIN data errors are computed using Gaussian error propagation producing an observational elevation change error (σ obs ). For the surfacefit and the crossover method the interpolation error is estimated as the RMS of the LSC uncertainty grid, defined as (σ int ). The final elevation change error is then estimated by combing the two error sources using RSS according to Here, Nρ is the number of uncorrelated grid cells estimated from empirical semi-variogram analysis of the CryoSat-2 and airborne elevation change differences, and estimated according to The Cryosphere, 10, 2953-2969, 2016 www.the-cryosphere.net/10/2953/2016/ where A mode is the total area covered by each CryoSat-2 mode and the correlation length ρ mode of 75 and 100 km for the LRM and SIN data, respectively.

Surface elevations compared to ATM
The measured surface elevations from the two CryoSat-2 products (JPL vs. ESA) showed large differences in both accuracy and precision of the elevation measurements, as seen in Table 1. The average accuracy and precision for the LRM data from the two products of 0.00 ± 0.43 m and −1.06 ± 0.89 m for the JPL and ESA products, respectively. This corresponds to an average reduction in RMSE of 68 % for the JPL product compared to the ESA LRM L2 data. Furthermore, our product shows a lower residual slope error (seen in Fig. 1c below ∼ 0.5 • ) indicating a lower sensitivity to the degradation of performance as the surface slope increases. Surface elevations generated from the SIN mode showed the same type of improvement as for the LRM data. Here, an average accuracy and precision were found to be −0.52 ± 0.58 m and −0.90 ± 1.05 m for the JPL and ESA SIN elevation products, respectively. This further corresponds to a reduction in the average RMSE of 27 % for the JPL product compared to the ESA product. For the SIN mode the JPL processing produces a slightly lower residual slope error, compared to the ESA processor (seen in Fig. 1c above ∼ 0.5 • ) Figure 2. Comparison of accuracy (a) and precision (b) as a function of the retracking threshold for the leading edge threshold retracker (LTH, dots) applied to the LRM (red) and SIN (blue) data and the leading edge maximum gradient retracker (LMG, dashed grey line) applied to the SIN data over Jakobshavn Isbrae and the region around the NEEM camp (77 • 27 N, 51 • 3.6 W) for the LRM data. The accuracy (mean) and the precision (standard deviation) has been determined for each threshold level using near-coincident ATM elevation, within a search radius of 50 m. The statistics were estimated using CryoSat-2 data from March to May 2013 and compared to ATM data from April 2013.
Larger improvements can be observed if separating the RMSE into its mean and standard deviation, corresponding to the accuracy and precision of the measurements. Using these definitions the analysis found that there is a 45 and 52 % increase in precision for the SIN and LRM data, respectively, compared to the ESA L2 product, and a 42 and 99 % improvement in accuracy for the respective modes.
The implementation of the LMG SIN retracking algorithm was found to reduce noise in the retrieved surface elevations compared to conventional threshold retracking. Though roughly comparable in accuracy, the LMG shows overall higher precision over all comparable leading edge thresholds. The adaptive nature of the algorithm provides improved estimates of surface elevation and a good trade-off between accuracy and precision.
The 20 % threshold retracker implemented in the LRM data was also found to provide improved estimates of surface elevation (both in accuracy and precision) compared to the model-based ESA L2 retracker. Furthermore, it also showed lower sensitivity to the 2012 melt event, due to the lower threshold used on the leading edge of the waveform.

Surface elevation changes compared to ATM
The estimated surface elevation changes generated from the surface-fit method also showed improvement in the estimated accuracy and precision, as seen in Table 2. Here, an overall improvement in RMSE of 55 and 40 % in the LRM and SIN data, respectively, was found when compared with ESA L2generated elevation changes from the same method. The average accuracy and precision, compared to ATM-generated elevation changes, were found to be 0.11 ± 0.67 m (LRM) and 0.30 ± 0.58 m (SIN) for the JPL-derived changes. This is compared to 0.25 ± 1.51 m (LRM) and 0.34 ± 1.06 m (SIN) www.the-cryosphere.net/10/2953/2016/ The Cryosphere, 10, 2953-2969, 2016 for the ESA-derived changes. This corresponds to an increase in elevation change accuracy of 56 and 12 % for the LRM and SIN data, respectively, for the JPL product compared to ESA L2 elevation changes. The estimated elevation changes for the JPL product also show an increase in precisions of 56 and 45 % for the LRM and SIN data, respectively, compared to its ESA counterpart. The estimated elevation changes of the Greenland Ice Sheet, excluding the peripheral glaciers, over the period of January 2011 to January 2015, show significant differences between products (JPL and ESA) in both spatial patterns and the total magnitude (Figs. 3, 4). The surface-fit and crossover methods produced on the order of ∼ 20 and ∼ 2.5 million usable elevation changes, respectively, providing high spatial sampling. Due to the constraint put into the JPL processor, the ESA L2 data produced slightly more surface-fit observations (∼ 10 %), as more surface elevations were accepted.
The ESA product produces a more positive elevation change pattern, which can be attributed to the 2012 melt event that introduced a large positive bias with a magnitude of ∼ 0.5 m (Nilsson et al., 2015a). Larger differences in the marginal areas for the surface-fit methods are also observed. The positive signal detected in the interior of the ESA surface-fit solution can also be found in the basin time series, correlating well with the timing of the summer of 2012 melt event, which can be seen in the time series in Fig. 3. These results are in agreement with earlier work demonstrating the sensitivity of the ESA retracker to the changes in the volume / surface scattering ratio (Nilsson et al., 2015a).
We used ICESat-and CryoSat-2-derived surface heights to generate time series over three regions in northeastern area of Greenland (Zachariae Isstrøm, Nioghalvfjerdsfjorden and Storstrømmen glaciers) for comparison purposes. These areas have recently shown large and rapid changes, which have been noted by, e.g., Khan et al. (2014). The selected areas were defined using hydrological basins derived by Lewis and Smith (2009), seen in Fig. 6, and were further divided into smaller areas around the termini to highlight performance for areas of rapid change. The ICESat and CryoSat-2 surface heights were then used to generate annual time series from 2003 to 2015 using Eq. (6). The estimated 12-year time series show overall comparable elevation change rates over both time periods (2003-2009 and 2010-2015), especially in the terminus areas, providing confidence that CryoSat-2 can reliably monitor changes in these areas.

Seasonal phase and amplitude compared to ICESat
The amplitude of the seasonal signal (Eq. 4) estimated from the surface-fit (SF) method shows large differences in both magnitude and spatial variability (Fig. 5). For the surface-fit method a difference in amplitude of 54 % is observed between the ESA and JPL products, corresponding to an areaaveraged amplitude of 0.17 m for the JPL product and of 0.37 m for the ESA product. The comparison with ICESatderived amplitudes from 2003 to 2009 estimated in Sasgen et al. (2012) using the same methodology as used here produced an area-averaged amplitude of 0.13 m, which is in good agreement with the JPL-derived amplitude. This agreement is also spatially consistent, as seen in Fig. 5, indicating low sensitivity to seasonal changes in the scattering regime of the upper snowpack. The observed difference in amplitude bias, taking ICESat as the true surface amplitude while acknowledging the differences in epochs and that no intercampaign bias has been applied, is 0.03 ± 0.13 m for the JPL product and 0. in amplitude is observed for all products. This is especially true for the ESA product, which increases its magnitude by a factor of 2. Analyzing the amplitude patterns on a regional drainage basin level (Fig. 5c) we find good agreement between JPL CryoSat-2 and ICESat amplitude with ESA data producing consistently larger amplitudes. Regionally, the highest amplitudes can be observed in southeastern Greenland in basins (3,4,5) and are consistent with regional precipitation patterns that show high total precipitation in these areas (Bales et al., 2009;Ettema et al., 2009).
The seasonal phase of the peak in amplitude of the seasonal cycle is shown in Fig. 5b, d and shows generally good agreement between the two data sets, providing the timing of the maximum of the accumulation signal, before the onset of melt, to the months of June/July for both JPL and ESA CryoSat-2 data sets. The ICESat-derived seasonal phase shows a higher dependence on elevation where the maximum of the accumulation signal is found in late May below 2000 m and late July/August above 2000 m in elevation. The ICESat discrepancies from the CryoSat-2 data are found in specific basins. Disagreements between the retrieved phase of the peak amplitude from CryoSat-2 and ICESat data are due to differences in temporal sampling as discussed in more detail in Sect. 8.

Volume change
The two volume change methods produce consistent results from JPL-derived elevation changes, with a difference of around 1 km 3 a −1 . The spread between volume change methods is larger (50 km 3 a −1 ) when using ESA L2 data. The larger discrepancy is mostly related to the sensitivity of the various methods to the melt event. The surface-fit method produces the most negative number (least affected by the melt event and has the lowest estimated error) and is therefore taken as the most reliable estimate for both the JPL and ESA solution.
Comparing the estimated volume change to other studies using CryoSat-2 we find that the JPL product is less negative than that estimated by Helm et al. (2014): www.the-cryosphere.net/10/2953/2016/ The Cryosphere, 10, 2953-2969, 2016 −375 ± 24 km 3 a −1 . This difference can be attributed to difference in processing methodology and to the different epoch of the data used by Helm et al. (2014) from January 2011 to January 2014. Using the corresponding epoch the JPL data gives a volume change estimate, based on the surface-fit method, of −353 ± 26 km 3 a −1 , well within the stated uncertainty of Helm et al. (2014).
To examine the regional behavior of volume change estimates of the Greenland Ice Sheet, gridded values from the two methods were divided into eight drainage basins according to Zwally et al. (2012). When analyzing the elevation time series at the basin scale, clear differences can be observed in the annual and inter-annual behaviors (Fig. 4). The northern and interior basins (1,2,7,8) all exhibit large differences (Table 4, 0-30 km 3 a −1 ) in the estimated volume change rates due to changes in the scattering regime resulting from the 2012 melt event. In the majority of the southern basins (4, 5, 6, 7), located in areas with higher precipitation, both products show good agreement in both trends and seasonal amplitude estimated from the surface-fit method.

Discussion
The CryoSat-2 processing methodology presented here is found to produce accurate and precise measurements of ice sheet elevation and elevation change. The main improvements have been introduced in the SIN processor with the inclusion of a novel type of land ice retracker (LMG), advanced phase filtering, and the inclusion of a phase ambiguity correc- tion scheme. This processing approach decreased the RMSE in the surface height retrieval by approximately 27 % (45 and 42 % improvement in precision and accuracy). This improvement further propagated into the quality of the estimated elevation changes for the SIN mode, with the same magnitude of improvement (Table 2). The described SIN processing also generated surface elevations and elevation changes with lower sensitivity to the local surface slope, indicating a higher degree of accuracy in the geolocation and surfacerange estimation. The SIN processing methodology further includes a phase filtering and phase ambiguity correction scheme. Visual inspections of a large number of tracks have shown more coherent estimation of the surface locations in the JPL product, and the implementation of the phase-ambiguity correction greatly reduced the number of track offsets. It was also noted that a relatively course DEM (∼ 1 km) could be used to resolve phase ambiguities. The detection and correction of phase ambiguities are relatively straightforward and rely mostly on the relative accuracy of the DEM. The implementation of the phase-ambiguity correction is particularly important when monitoring smaller ice caps and outlet glaciers, where frequent and large track offsets can bias the estimation of the underlying topography.
The new LRM data processing methodology focused on optimal retrieval of surface elevations over the interior parts of the ice sheet. Here the choice of retracking threshold has proven to be the critical factor to acquire high-quality surface elevations and elevation changes. The choice of 20 % leading edge threshold level reduced the sensitivity to changes in the scattering regime for low-slope, high-elevation areas. The functional-based retracking algorithm used in the ESA LRM data processor corresponds roughly to a 50 % threshold level (Wingham et al., 2006a), which appears to suffer from a higher sensitivity to changes in the scattering properties (volume scattering) of the near-surface firn, as the range is referenced higher up (later in time) on the leading edge of the waveform. This effect can be seen in Fig. 2a. It can also The Cryosphere, 10, [2953][2954][2955][2956][2957][2958][2959][2960][2961][2962][2963][2964][2965][2966][2967][2968][2969]2016 www.the-cryosphere.net/10/2953/2016/ be seen that the observed negative elevation bias (Table 1) for ESA LRM (−1.0 m) fits well with the bias for the 50 % LRM data threshold value shown in Fig. 2a. This makes the algorithm more sensitive to annual and sub-annual changes in snowpack volume / surface scattering ratio, which can produce spurious changes in elevation due to changes in the near-surface dielectric properties. This is clearly shown in patterns of ESA-product-derived elevation changes (Fig. 3) where a large elevation bias was introduced by the 2012 melt event (Nilsson et al., 2015a). The 20 % threshold is less sensitive to these types of changes (Tables 1, 2) and the results are in agreement with previous work that has demonstrated that the 20 % threshold best represents the mean surface inside the footprint when exposed to a combination of surface and volume scattering (Davis, 1997). Surface elevation changes, derived from multi-temporal radar altimetry observations, are typically corrected for their correlation to changes in the radar waveform shape. This is to reduce the effect of changes in the volume / surface scattering ratio of the ice sheets' surface (Davis, 2005;Flament and Rémy, 2012;Wingham et al., 2006b;Zwally et al., 2005). This inherently adds to the complexity of the processing and analysis, introducing new biases and error sources in the estimated parameters. For the processing approach presented here many of the post-processing steps can be omitted or reduced, as they are an inherent part of the improved waveform retracking. There have been attempts to remove spurious step changes in elevation resulting from sudden changes in surface-scattering characteristics (caused by the 2012 melt event) apparent in the ESA Baseline-B L2 data through postprocessing strategies (Nilsson et al., 2015c;McMillan et al., 2016), but such approaches spread the bias over a longer period of time making the "jumps" less noticeable in the time series by removing the step change. However, these strategies introduce longer-timescale bias of equal magnitude as the scattering layer is buried by less reflective snow and lowdensity firn.
The result of the validation procedure shows a larger slope-dependent bias in the ESA data, both in the elevation and elevation changes (Fig. 1). This is especially true for the surface elevations, which can be seen in the figures of precision and accuracy (Fig. 1a, c), where both figures show clear linear slope for the ESA surface heights. In comparison, estimated elevations from JPL product show relatively stable statistics over the entire slope range above 0.2 • . The validation of the estimated surface elevation changes, seen in Fig. 1b and d, shows the effect of the 2012 melt event on the ESA-derived elevation changes below 0.2 • . Furthermore, the accuracy of the ESA-derived changes shows a clear negative trend as function of increased surface slope. The derived precision of the surface elevation change decreases dramatically above 0.5 • , as more complex topography is measured.
The JPL CryoSat-2 processing methodology produces seasonal amplitudes that are in good agreement with those derived from ICESat data, further indicating the processors' abilities to track real and physical changes of the ice sheet surface. The current ESA implementation produces noisier estimates of elevation change, as indicated by the larger standard deviations of the residuals in the ESA solutions for the surface-fit and crossover methods. Figure 5 further shows an amplitude bias in the ESA data compared to the corresponding ICESat reference amplitudes. The bias is constant above the Greenland ELA, located around 2000 m in altitude, but increases linearly as elevations decrease below this. The linear increase in amplitude seems to be connected to the higher and more variable precipitation in the ablation zone where changes in the variable snow cover produce changes in apparent surface height. This is less prominent for the JPL SIN and LRM data retrackers. The estimated seasonal phase in Fig. 5b and d shows that both JPL and ESA CryoSat-2 elevation products can adequately resolve the seasonal maximum of the accumulation signal. Both products provide a timing of the maximum to the month of July over the entire ice sheet, independent of elevation. Assessing the CryoSat-2-derived maximum reveals a difference between CryoSat-2 and the reference ICESat data set. This constitutes roughly a ±1-month difference depending on the elevation and the location. The cause of this difference can be attributed to the temporal sampling of the ICESat mission. During the mission, due to degraded laser lifespan, data were only collected in campaign mode during the spring and winter times corresponding to roughly 2 months of measurements for each period. When the CryoSat-2 data were resampled to coincide with the ICESat temporal sampling the same elevation and spatial pattern in the phase of the maximum seasonal amplitude was observed as determined from the ICESat data. No corresponding change in amplitude was observed. This was done by selecting CryoSat-2 data corresponding to the same unique months available in the total ICESat record.
The two independent methods used to estimate the volume change of the Greenland Ice Sheet produce consistent volume change estimates. This was especially true for volume changes derived from the JPL elevations, with a discrepancy of less than 1 km 3 a −1 between methods. The two methods provided the same estimate of integrated volume change, but the use of the surface fit is recommended as it produces higher spatial sampling compared to the crossovermethod and lower errors. The good agreement between the methods further indicates a strong reliability in the estimated volume change rates of the Greenland Ice Sheet over the 4-year period. It also shows the ability of CryoSat-2 to capture both small-and large-scale spatial patterns in the rugged topography along the coastline and in the interior of Greenland. This is especially true in the major outlet-glacier systems (e.g., Zachariae Isstrøm, Nioghalvfjerdsfjorden, and Storstrømmen).
Studying the northern parts of the Greenland Ice Sheet we find that CryoSat-2 captures both intricate and complex behavior in the marginal areas of the ice sheet. This is exemplified in the northeastern regions of Greenland (Fig. 6)  Zacahariae Isstrøm, Nioghalvfjerdsfjorden, and Storstrømmen, which all show complex and localized patterns of elevation change. Here, Nioghalvfjerdsbrae shows very small changes in elevation during the observational time span, while Zacahariae Isstrøm, its major neighbor, shows large negative trends in elevation change. The observed behavior agrees with the observations made in recent studies by Khan et al. (2014) and Mouginot et al. (2015), who document rapid retreat and drawdown of the ice-front position of the Zacahariae Isstrøm beginning in 2012. Storstrømmen outlet glacier system also appears to show signs of rapid thinning at low elevations near the ice-front position, while a large positive signal is observed roughly 100 km upstream of the terminus. This pattern has also been observed by Joughin et al. (2010) and Thomas et al. (2009), using airborne altimetry and surface velocity mapping. Rates of elevation change from ICESat and CryoSat-2 data show good agreement in basin-scale trends (Fig. 6b, c)   . Converting the estimated mass change to volume change, assuming no changes in firn air content over the study period and an ice density of 917 kg m −3 (assessment of changes in firn air content is out of the scope of this paper), gives an ice sheet wide rate of −305 ± 38 km 3 a −1 (updated to Schlegel et al., 2016). This estimate is corrected for volume changes from peripheral glaciers that lost volume at a rate of approximately −40 ± 27 km 3 a −1 (Box, 2013;Fettweis et al., 2013;Noël et al., 2015). This estimate of peripheral glacier change is in agreement with the estimated volume change of −41 ± 8 km 3 a −1 from Gardner et al. (2013). The volume rate derived from GRACE data agrees well with our estimated rate from CryoSat-2, where both results are within the 1σ uncertainty of each other and neglecting changes in firn air content over the period of study. The observed volume change rates estimated from this study are within the range of previous studies, ranging from −186 to −309 km 3 a −1 for the time period 2003-2009, summarized by Csatho et al. (2014), using the same mass-to-volume conversion applied to the GRACE data. A more recent study by Helm et al. (2014) of −375 ± 24 km 3 a −1 agrees within uncertainties when differences in observation periods (2011-2014 vs. 2011-2015) are taken into account. From this comparison we find that our estimate spans both the estimate of Csatho et al. (2014) and the mass loss estimated from GRACE, while acknowledging the varied time spans of the different studies.

Summary and conclusion
We conclude that the use of an adaptive retracker for the SIN mode, based on the maximum gradient method, and the use  (2011-2015), with corresponding hydrological basin outlines. The hydrological basins are separated into full basin size (b) and to the terminus areas (c). (b, c) Show a merged 12-year annual elevation time series from ICESat and CryoSat-2 for each color-coded area in (a). The derived elevation time series was formed using the surface-fit method described in Sect. 3.1 onto a 500 m grid to facilitate merging of the two data sets due to their difference in orbit characteristics. The elevation change map is overlaid onto the CryoSat-2 hillshaded DEM based on surface heights from July 2010 to February 2015. The annual 12-year time series was created from the surfacefit method by binning the monthly values into annual values using the median of the corresponding 12 months. The gray vertical line indicates the temporal separation between the ICESat and CryoSat-2 missions. of a 20 % threshold retracker for the LRM data provide improved performance to the retracker currently used for the ESA L2 elevation products. It is further important, especially for the SIN mode, to apply a leading edge discriminator to identify and track the leading edge of the waveform. The functional model currently employed in the ESA processor has, to the author's knowledge, no such discriminator currently implemented. This is important in the SIN mode, as it often contains multiple surface returns. The single-return model applied in the ESA processor will have issues fitting a waveform containing multiple surface returns resulting in retrack jitter (Helm et al., 2014).
Using the new CryoSat-2 processing methodology for the LRM and SIN data we determine the volume change of the Greenland Ice Sheet to be −289 ± 20 km 3 a −1 during the period January 2011 to January 2015. The validation against airborne ATM surface elevations and elevation changes showed an average improvement in the RMSE of the measured elevations of 68 and 27 % for the LRM and SIN data, respectively, compared to ESA Baseline-B L2 products. The new methodology also provides improved elevation The Cryosphere, 10, 2953-2969, 2016 www.the-cryosphere.net/10/2953/2016/ changes with an reduction in RMSE of 55 and 40 % for the LRM and SIN data, respectively, compared to their ESA L2derived data.
The methodology also showed less sensitivity to changes in near-surface scattering properties than equivalent ESA products. The new processing methodology showed little effect of slope-induced errors, providing better performance in the marginal areas of the ice sheets. These improvements to the CryoSat-2 processing mitigate the need for postprocessing to correct for the correlation between changes in surface elevation and changes in the waveform shape (i.e., backscatter and leading edge width etc.) that can introduce biases and add to the complexity of the processing and analysis.
The presented CryoSat-2 processing methodology provides a lower intrinsic error in the measured elevation, elevation change, and volume change estimates, all of which will facilitate improved understanding of the geophysical process leading to changes in land ice elevation. Given the release of the ESA Baseline-C, which provides improved corrections and processing mainly for the L1B product, further improvements are expected in the near future.

Data availability
The complete DEM and elevation change grids used in this study are available upon request to the corresponding author and are provided in GeoTIFF format The Supplement related to this article is available online at doi:10.5194/tc-10-2953-2016-supplement.