An elevation change dataset in Greenland ice sheet from 2003 to 2020 using satellite altimetry data

ABSTRACT A decade-long pronounced increase in temperatures in the Arctic resulted in a global warming hotspot over the Greenland ice sheet (GrIS). Associated changes in the cryosphere were the consequence and led to a demand for monitoring glacier changes, which are one of the major parameters to analyze the responses of the GrIS to climate change. Long-term altimetry data (e.g. ICESat, CryoSat-2, and ICESat-2) can provide elevation changes over different periods, and many methods have been developed for altimetry alone to obtain elevation changes. In this work, we provided the long-term elevation change rate data of the GrIS in three different periods using ICESat data (from February 2003 to October 2009), Cryosat-2 data (from August 2010 to October 2018) and ICESat-2 data (from October 2018 to December 2020). Optimal methods were applied to the datasets collected by three different altimeters: crossover analysis for ICESat/ICESat-2 and the surface fit method for Cryosat-2. The data revealed that the elevation change rates of the GrIS were –12.19 ± 3.81 cm/yr, –19.70 ± 3.61 cm/yr and –23.39 ± 3.06 cm/yr in the three different periods, corresponding to volume change rates of –210.20 ± 25.34 km3/yr, –339.11 ± 24.01 km3/yr and –363.33 ± 20.37 km3/yr, respectively. In general, the obtained results agree with the trends discovered by other studies that were also derived from satellite altimetry data. This dataset provides the basic data for research into the impact of climate change over the GrIS. The dataset is available at https://doi.org/10.57760/sciencedb.j00076.00121.


Introduction
Over recent decades, the Greenland ice sheet (GrIS), which is the second-largest ice sheet in the world, has experienced substantial mass loss driven by both discharge and surface melt (Hanna, Mernild, Cappelen, & Steffen, 2012;Straneo & Heimbach, 2013).The Sixth Assessment Report of the Intergovernmental Panel on Climate Change pointed out that the mass loss of the GrIS reached 4,890 (4,640) Gt from 1992 to 2020, equivalent to 13.5 (11.4-15.6)mm global mean sea-level rise, which has led to the bedrock of Greenland rising tens of centimeters from 2007 to 2019 (Gulev, Thorne, & Vose, 2021).The surface elevation change in the GrIS is a direct indicator of climate change and ice dynamics.Additionally, it can be used to derive mass changes of the GrIS with the consideration of the surface mass balance, firn compaction, dynamics and underlying bedrock motion, and can validate prognostic models simulating the recent evolution of ice sheets (Shepherd et al., 2012;Zwally et al., 2011).Therefore, long-term and large-scale GrIS elevation change detection is of great significance for the study of global climate change.
Satellite altimetry is an effective technology to detect changes in the GrIS, which can compensate for the limitations of spatial and temporal coverage of airborne and field observations.It transmits the laser or radar signal with known power toward the nadir direction and measures the traveling time of the signal from the altimeter to the ground surface and back.The elevation observation can be estimated by the time and free-space speed of light (Chelton, Ries, Haines, Fu, & Callahan, 2001).Since the early 1990s, a series of satellite radar altimetry has been used to measure elevation changes in Greenland such as ERS-1, ERS-2, Envisat and currently Cryosat-2 (Arthern, Wingham, & Ridout, 2001;Davis, Kluever, & Haines, 1998;Sørensen et al., 2015;Zwally, Bindschadler, Brenner, Major, & Marsh, 1989).Radar signals with longer wavelengths are insensitive to cloud cover and can therefore enable continuous observations.For laser altimetry, NASA's ICESat in a period between 2003 and 2009 and ICESat-2 from 2018 are the two main data sources for the estimation of ice sheet elevation.Compared to radar altimeters with up to a kilometer footprint radius on a flat water surface, laser altimetry has a smaller footprint and lower slope-induced errors as well as repeat track accuracy (Brenner, DiMarzio, & Zwally, 2007;Fricker, 2005).The single-point accuracy of the measured elevation is usually within 0.1 m to 0.15 m, and the distance between repeat tracks is within a few hundred meters (Shuman et al., 2006).However, laser measurements are affected by clouds or drifting snow because of its short wavelength, and the ICESat spans only six years (Arthern, Wingham, & Ridout, 2001;Ridley & Partington, 1988).ICESat-2, as a continuation and improvement of ICESat, is providing an invaluable dataset of global elevation measurements to support diverse scientific applications due to its dense spatial and temporal coverage.In addition, CryoSat-2 is a suitable radar satellite altimeter system providing surface elevation data of the GrIS from 2010 to 2018, when there was a gap in observations between the ICESat and ICESat-2 missions.Therefore, the combined application of ICESat, CryoSat-2 and ICESat-2 provides multiperiod elevation change information for the GrIS.
Multiple methods have been developed for different altimetric data to retrieve the elevation change.The most common ones are crossover, repeat track and plane fit analysis.In some early studies, Johannessen, Khvorostovsky, Miles, and Bobylev (2005) used crossover analysis to estimate the elevation change rate from 1992 to 2003 based on ERS-1 and ERS-2 data.Khvorostosvsky (Ewert, Groh, & Dietrich, 2012) added Envisat data into this method and extended the timespan to 2008.In addition, Ewert, Groh, and Dietrich (2012) developed the repeat track method for ICESat data to calculate the elevation change, volume change and mass balance of the GrIS from 2003 to 2008.In recent studies, Sørensen et al. (2018) applied a combination of crossover, repeat track and plane fit algorithms to ERS-1, ERS-2, Envisat and Cryosat-2 radar data and presented 25 years of elevation change estimates (1992-2016) of the GrIS.Levinsen et al. (2015) compared the crossover and along track results over the Jakobshavn drainage basin and concluded that the crossover method was advantageous for ignoring slope-induced errors and that the along-track method was advantageous for large spatial data coverage.In this study, we produced an elevation change dataset using ICESat data from February 2003 to October 2009, CryoSat-2 data from August 2010 to October 2018 and ICESat-2 data from October 2018 to December 2020, which are widely accepted due to the higher accuracy of the three data (Brenner, DiMarzio, & Zwally, 2007;Pritchard, Arthern, Vaughan, & Edwards, 2009;Shuman et al., 2006).We updated the data span time of obtaining elevation changes based on previous studies and analyzed the elevation and volume changes of GrIS over a longer period of time.We also estimated the elevation change rate during the operational period of the ICESat-2, which was less in the current studies.To avoid system error between missions, we present the elevation change results from the three data separately.
The dataset is available at http://doi.org/10.57760/sciencedb.j00076.00121.The results were analyzed using a drainage scale and compared with the volume loss from previous studies.This will allow the dataset to be further utilized in other research studies over the GrIS and integrated with multiple datasets.

Icesat and ICESat-2 data
The ICESat satellite, the first laser altimeter satellite mission for Earth observation, was launched by NASA on January 12, 2003 and operated 18 measurement campaigns from 2003 to 2009.It hosted the Geoscience Laser Altimeter System (GLAS) payload aiming to measure ice sheet mass balance and cloud and aerosol heights, as well as land topography and vegetation characteristics.These data demonstrated a remarkable capability to assess changes in ice sheet elevations despite its compromised operations (Schutz, Zwally, Shuman, Hancock, & DiMarzio, 2005).
In this study, we used the GLA12 release 34 products provided by NSIDC from February 2003 to October 2009 to estimate the elevation change in the GrIS.We carried out the following criteria data-culling and correction procedures to reduce systematic errors and outliers: 1) sorted out data not meeting the on-product variable criteria: the quality flag is valid (elev_use_flg = 0), saturation correction is applicable (sat_corr_flg < = 2), attitude quality indicates good (sigma_att_flag = 0) and the number of peaks in the waveform produced is 1 (i_numPk = 1); 2) applied the inter-laser correction: subtracting 1.7 cm from Laser 2 and adding 1.1 cm to Laser 3 (Smith et al., 2020); 3) converted the elevation of the Topex/Poseidon ellipsoidal to elevation of the WGS84 ellipsoidal to align with the coordinates of ICESat-2 and Cryosat-2.
As the successor of ICESat, ICESat-2 was launched on 15 September 2018 with the Advanced Topographic Laser Altimetry System (ATLAS) instrument to measure ice sheet mass changes, cloud and aerosol heights, land topography, and vegetation characteristics (Pattyn et al., 2018).To measure the Earth's surface height more accurately, the ICESat-2 laser is split into six individual beams in three pairs.The beams in a pair have different transmit energies (energy ratio of ~1:4) and are separated by 90 m in the across-track direction.The key mission parameters of ICESat and ICESat-2 are shown in Table 1.
Compared with the laser on ICESat that sends 40 pulses per second, the ICESat-2 laser is fast-firing and sends 10,000 pulses per second.ICESat-2 has denser footprints than ICESat.
The multibeam instrument design, smaller footprint, and the ability to resolve rougher terrains enable more accurate elevation measurements and fewer measurement gaps (Abdalati et al., 2010;Markus et al., 2017).The ICESat-2 Land Ice Height ATL06 dataset, which was developed from ATL03 photon data from October 2018 to December 2020, was used for our dataset (Smith et al., 2019).This includes latitude, longitude, and height, with respect to the WGS84 ellipsoid for each segment.Only data marked as good quality (atl06_quality_summary = 0) were used.

CryoSat-2 data
The CryoSat-2 radar altimetry satellite, the first ice mission of ESA, was launched on April 8, 2010.The main payload is an advanced radar altimeter, the SAR interferometric radar altimeter (SIRAL), which was designed to detect changes in sea ice thickness as well as surface elevation changes in Earth's land and marine ice sheets.It operates in two modes over the GrIS: low-resolution mode (LRM), which is operated over the flat interior part, and SAR interferometric mode (SARIn), which is operated over the complex steep terrains on the edges.Its key parameters are shown in Table 1.Its narrow across-track spacing of 2.5 km to 4 km at the margins of the ice sheet is an improvement over the ICESat instrument's 25 to 40 km spacing.
In our dataset, we employed the Level-2 LRM and SARIn Baseline-D product obtained from August 2010 to October 2018 to estimate the elevation change of GrIS.

Data for validation
The Operation IceBridge Airborne Topographic Mapper (ATM) LiDAR system has operated every spring in Greenland since 1993.The ATM data provide a valid comparison for satellite data because laser altimeters generally have better accuracy than radar altimeters.In addition to the surface elevation (Level 2) product from multiple campaigns, a Level 4 (L4) product is provided which consists of surface elevation changes derived from repeated observations in the ATM dataset.We used ATM L4 from 2003 to 2009 and ATM L4 from 2010 to 2018 as an independent validation dataset to evaluate our elevation change results obtained from ICESat and CryoSat-2.A 40 km floating median low-pass filter was used to eliminate outliers in the ATM L4 data before performing the validation (Zhang, Wang, An, Liu, & Geng, 2022).

Elevation variation estimation using ICESat and ICESat-2
For ICESat and ICESat-2 data, the crossover analysis method was applied.Because of the design characteristics of the ICESat/ICESat-2 satellite orbit, ascending and descending tracks from different campaigns always intersect.The crossover method was used to derive the surface elevation change at the crossover between an ascending and descending satellite ground track separated in time (Abdalati et al., 2010;Brenner, DiMarzio, & Zwally, 2007;Smith, Bentley, & Raymond, 2005).
The elevation and time at the crossover can be linearly interpolated using the two closest observations along the track.Two interpolated elevation values from two different interpolated epochs for each track crossover are expressed by Eq. ( 1).The interpolated elevations are subtracted, yielding an elevation difference (H), and the interpolated measurement times are subtracted, yielding a time difference (t).
where H a , H b , H c and H d are observations along the different tracks.A 1 , A 2 , D 1 and D 2 are the corresponding distances to the crossover.
To reduce the elevation error caused by the slope and the interpolation distance, the observations used for interpolation should not exceed 550 m for ICESat and 10 m for ICESat-2 data (i.e., no more than two missing observations along the track).This method can reduce errors due to terrain slope, and shows higher precision (Felikson et al., 2017).Thus, only track points that meet the requirements were used.Figure 1(a) shows the crossover distribution of ICESat, and Figure 1(c) shows the number of crossover for ICESat-2 data in each 5 km grid.The crossover points are denser in high latitude areas but less dense in low latitude areas, whereas ICESat-2 generated 20 times as many crossover points in two years as it did in the entire ICESat mission.
Due to the errors in the satellite control system, the tracks of the laser altimeter are not precisely repeated and the crossover of ground tracks can be up to hundred meters apart.Assuming that ice sheet elevation changes are consistent within a 5 km grid area, these intersections can be used to determine the rate of elevation change.The local ice surface elevation within each grid cell can be expressed by where H t ð Þis the elevation at crossover.H 0 x; y ð Þ is the topographic term and H 1 t ð Þ is the fluctuation term.dh dt is the long-term linear ice elevation change rate.asin ωt ð Þ and bcos ωt ð Þrepresent the periodic fluctuation in the annual elevation difference.According to Eq. ( 2), the elevation difference at crossover can be expressed by Eq. ( 3), where the topographic term H 0 x; y ð Þ can be offset.
where H is the elevation difference between the interpolated elevation along the ascending and descending tracks, and t 1 andt 2 are the corresponding interpolated times.According to Eq. ( 3), we performed a least-squares regression to calculate three parameters, dh dt ,a and b within each grid cell.To reduce the effect of any poor fit, quality control was constrained in terms of data rationality.We computed the standard deviation of residuals and excluded individual measurements whose corresponding residuals exceeded three times of the standard deviation (3σ).We iteratively repeated the parameter fit until no more outliers were found.Furthermore, to exclude the remaining unrealistic results from further processing, we rejected any results where we obtained an absolute elevation change rate ( dh dt ) larger than 10 m/yr or where the standard deviation of this rate was higher than 0.5 m/yr (Khan et al., 2014;Schröder et al., 2019;Slater et al., 2018).We used a minimum time span of 2 years for ICESat and 1 year for ICESat-2 to minimize short-term meteorological variations in the dh dt (Sørensen et al., 2011).

Elevation variation estimation using CryoSat-2
The orbit of Cryosat-2 is not necessarily located in a straight and repeat line and may be more scattered than that of the conventional laser altimeters.The surface fit method has been proven to be a better way for CryoSat-2 data in existing studies to determine surface elevation change (Nilsson, Gardner, Sandberg Sørensen, & Forsberg, 2016).In our work, this method was used.
Surface fit was performed by fitting a least-squares model to the elevations as a function of time and space inside a relatively small area under the assumption that elevation changes were consistent in a 5 km grid area (Simonsen & Sørensen, 2017).Only observations from the ice sheet were needed, and the number of observations for each 5 km grid is shown in Figure 1(b).As mentioned above, the elevation at a grid can be expressed as a combination of the topographic and fluctuation terms.For Cryosat-2, the topographic term (H 0 ð� x; � yÞ) was fitted by the cubic polynomial.In addition, a backscatter correction factor and the leading edge of the radar echo were applied to alleviate the impact of changes in ice sheet surface penetration.Furthermore, the bias between ascending and descending (AD) orbits should not be ignored (Legresy, Frédérique, & Philippe, 1999;Armitage et al., 2013;McMillan et al., 2014).Hence, the surface fit can be expressed as follows: where H (x, y, t) is the elevation measured by Cryosat-2.H 0 is the mean altitude in the grid cell.c 1 to c 4 are local topography coefficients.dh dt is the elevation change constant and the superposition of sine and cosine represents the periodic fluctuation of annual elevation difference.B s and LeW are the backscatter coefficient and the leading edge width, respectively.b AD is the ascending and descending bias, and AD is the logical index (0/1) describing the flight direction of the satellite.b m is the CryoSat-2 mode bias and m is the logical index (0/1) describing the model.
For CryoSat-2, we applied the same 3σ criterion and threshold in the iterative calculation as the least square fitting of ICESat/ICESat-2 (Eq.( 3)) for quality control.In addition, we rejected the solutions if the b m was greater than 1 m or the seasonal variation in elevation ( ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ) was estimated to be greater than 3 m (Felikson et al., 2017).Furthermore, we assumed that the root mean square error (RMSE) of the modeled elevation was 10 m (Fan, Ke, & Shen, 2022).The minimum number of observations and the minimum timestamp were set to 30 and two years respectively for Cryosat-2.

Uncertainty quantification
For the calculated grids, least-squares adjustment produces not only unique values for the unknown parameters but also indicates their precision.The regression process for elevation change can be written as a matrix equation: where Y is a vector of observations.B is the design matrix of Eq. (3) for ICESat/ICESat-2 and Eq.(4) for CryoSat-2.X is the parameter vector to be estimated, and V is the vector of residuals.
To illustrate the uncertainty of the estimated change rate in this study, we adopted the variance -covariance matrix (VC-matrix) Q xx as the error evaluation indicator (Slobbe, Lindenbergh, & Ditmar, 2008), which is expressed in Eq. ( 6).
where B is the design matrix relating the measurements to be estimated and P is the weight matrix of the height difference, which we assume to be 1.S 0 is the unit weight standard deviation.n is the number of independent observations in the grid cell and t is the number of parameters to be solved.Q dh/dt is the square root of the diagonal element of the variance -covariance matrix Q xx .S dh/dt is the uncertainty of the elevation change rate.
For the uncertainty estimation of interpolated grids, we used the kriging variance error calculated with ArcGIS 10.6 (Fan, Ke, & Shen, 2022).There is a 95.5% probability that the actual elevation at the grid is the predicted raster value ±2 times the square root of the variance error of the corresponding cell by assuming that the kriging errors are normally distributed.We take two times the square root of the value in the variance error as the elevation change uncertainty in the interpolated grids (Eq.( 7)).
where e is the kriging variance error.

Elevation change and uncertainty during three periods
Due to the different time spans of altimeter data, we produced the GrIS elevation change rate dataset during three different periods using ICESat, Cryosat-2 and ICESat-2 data.The interpolated 5 km × 5 km grids from the surface elevation change rates using ordinary kriging are given in Figure 2, with red regions in areas of thinning while blue regions are thickening.We applied a 5 × 5 median filter over all regions to remove outliers in the gridded dh dt values.The spatial pattern of the elevation change rate during the three periods from ICESat, CryoSat-2 and ICESat-2 showed a good agreement.A distinct thinning signal of approximately 1.5 m/yr could be found along the GrIS margin, and a relatively lower thinning or stable thickening was found in the interior part of the ice sheet.The average elevation loss rates derived from ICESat, Cryosat-2 and ICESat-2 were -12.19 ± 3.81 cm/yr, -19.70 ± 3.61 cm/yr and -23.39 ± 3.06 cm/yr, respectively.The highest thinning rates occurred at the Jakobshavn and Kangerlussuaq glaciers, with rates of 4 and 6 m/yr, respectively, while the largest inland thickening rate was 0.15 m/yr.We derived the volume changes by multiplying the grid cell size (5 km) with the corresponding corrected elevation changes.In general, the volume loss during the Cryosat-2 period from 2010 to 2018 was -339.11km 3 / yr, approximately 128.91 km 3 /yr more than that during the ICESat period from 2003 to 2009 (-210.20 km 3 /yr).The loss during the ICESat-2 period increased by approximately 24.22 km 3 /yr compared to the Cryosat-2 period.The negative values indicate that the GrIS was still in a state of obvious mass loss as a whole from 2003 to 2020.
In Section 3.3, we described the approach to estimating the uncertainty in our results.Figure 3 shows the uncertainty of the estimated elevation change rate from ICESat, Cryosat-2 and ICESat-2 data in each grid.The uncertainty of elevation change was larger in the marginal region of the GrIS due to the influence of topography, which could be as large as 0.50 m/yr over some areas.In contrast, the inland area was smaller due to more observations and flat terrain, which was approximately 0.01 m/yr for CryoSat-2 and ICESat-2 data (Smith et al., 2020).However, the interpolated or extrapolation process to obtain the elevation of the crossover point can introduce errors in interpolated grid cells.This error is more obvious for the ICESat data.As the ICESat result have a larger data gap than ICESat-2 and Cryosat-2 results due to sparse observations, especially at low latitudes, the interpolated error can range from approximately 0.04 to 0.05 m/yr.

Elevation and volume change rates in different drainage basins
The elevation and volume change rates during the three periods in each drainage basin are given in Table 2.In addition, Figure 4 shows the volume change rates of the GrIS drainage basin divided by the elevation contour line of 2000 m.
The SE drainage had an overall severe and accelerated volume loss from ~-100 km 3 /yr to ~-140 km 3 /yr, in addition to a significant thinning rate up to ~-4 m/yr at elevations <2000 m.The NW drainage exhibited a much smaller volume decrease of ~-60 km 3 /yr at elevations <2000 m.The highest rates of decrease could be found around the marine terminal west of this drainage.In the northern GrIS, the NO drainage experienced the smallest overall thinning rate, ~-0.1 m/yr, and a volume decrease of ~-30 km 3 /yr at altitudes under 2000 m.Thickening at higher altitudes at a much smaller magnitude was observed (Hanna et al., 2021).In the SW drainage, measured thinning of up to ~-40 km 3 /yr was found at elevations <2000 m.CW was undergoing stable ablation of approximately -40 km 3 /yr.The glaciers in CW with the accelerated speed of elevation decline indicated more serious deglaciation (Howat, Smith,  Joughin , & Scambos, 2008).Hotspots of decrease were found at the Jakobshavn glacier up of ~-4 m/yr.In contrast, the NE drainage was the only region that had significant accumulation, ~40 km 3 /yr at elevations >2000 m.Slight thickening at lower altitudes could be found, while thinning was limited to the margins, especially in the Zachariae glacier, up to ~-3 m/yr.

Comparison with airborne laser altimetry
In this section, the elevation change results of the ICESat and CryoSat-2 data are compared with the Operation IceBridge ATM L4 elevation change product.The ATM elevation changes, the difference in the elevation change rate between ICESat/Cryosat-2 and ATM are shown in Figure 5.The histogram of elevation change differences at 0.1 m/yr bins is in the lower left corner of Figure 5(b,d).In each grid cell, the median of all IceBridge data was subtracted from our results.Figure 5(b,d) show the elevation change difference, the red and blue spots on the outlet glaciers indicate a serious overestimation and underestimation, respectively.We can see that the larger differences are mainly concentrated on the edge of the GrIS where is with complex terrain.Especially over the Jakobshavn glacier basin, our results slightly overestimate the elevation change rate.Subsequently, the median and root mean square (RMS) statistics of the differences for different specific surface slopes are summarized in Table 3.The surface slope was obtained from an ICESat-2-derived digital elevation model (DEM) (Fan et al., 2022).The median difference in elevation changes between ICESat and ATM amounted to -0.12 m/yr with  The IceBridge data cover the margins of many GrIS glaciers, where elevation changes differ over relatively short distances.We found a significantly larger spread of RMS in higher slope areas.In flat areas with the surface slope <0.2°, the two results have the difference of only ~0.02 ± 0.4 m/yr with ATM L4.With the increase in surface slope, the differences with ATM increase rapidly and reach ~0.30 ± 1.80 m/yr when the surface slope is more than 1°.Possible reasons for this may include measurement errors, poor observation accuracy for satellite altimeters in the complex terrain areas, and interpolation errors (Lai & Wang, 2021).In the flat regions, the RMS of the differences is smaller, at a level of 10 cm/yr.

Comparison with previous studies
To further evaluate the performance of our data, Table 4 shows the comparison of GrIS volume change rates between this study and the other four published studies (Felikson et al., 2017;Simonsen & Sørensen, 2017;Smith et al., 2020;Li et al., 2022).In general, our estimations show an overall good agreement with their results.
For ICESat, the GrIS volume change rates over the same period calculated by this study were higher than that in Felikson et al. (2017).Nonetheless, a diffidence of ±50 km 3 /yr is within the uncertainty of glacier volume change.This is potentially attributed to the different regression models used.Although we both used crossover analysis, the leastsquares regression used by Felikson et al. (2017) was a linear model as follows: where it contains a linear elevation change rate and the error term.Our study for crossover analysis contains an extra periodic fluctuation of the annual elevation difference term as shown in Eq. (3).
For CryoSat-2, we used the same method as LSM4 in Simonsen and Sørensen (2017).Their study obtained a volume loss of -271 ± 32 km 3 /yr between 2010 and 2014 for the GrIS.We updated the span time to 2018 for CryoSat-2, and our values show that there was an increasing volume loss of the GrIS of approximately -30 km 3 /yr.In more recent studies, Smith et al. (2020) and Li et al. (2022) estimated trends of -235 ± 3 km 3 /yr and -192 ± 7 km 3 /yr, respectively, from 2003 to 2019 utilizing ICESat and ICESat-2.They generated the crossover between ICESat and ICESat-2 and estimated the elevation surface changes over a longer period of time.So that they obtained a less volume loss rate.

Limitations of our dataset
We used a series of processes to ensure the reliability of our results and showed an empirical estimation of the uncertainty distribution in Section 4.1.However, according to previous studies, there still exist some factors that might introduce errors into our results.These factors are important but difficult to formally quantify.First, variation in penetration depth into the snow is a typical characteristic of the radar altimeter, which causes Cryosat-2 to not observe the actual surface elevation change of the ice sheet.Previous studies have shown that the penetration depth in a dry snow zone can reach several meters (Slater et al., 2019).Under stable climate conditions, this might actually be an advantage, as it might suppress noise induced by seasonal snowfall and give an elevation change estimate more relevant than that from laser altimetry (Sørensen et al., 2018).However, the past decade has been characterized by large climate variability in Greenland, and the surface penetration is highly impressionable.For example, the elevation record of Cryosat-2 during the strong melt event in July 2021 showed that an abrupt increase resulting from surface melting created strong reflective surfaces for the radar signal (Forsberg, Sørensen, Levinsen, & Nilsson, 2013;McMillan et al., 2016;Nilsson et al., 2015).Therefore, some surface processes, such as melting and refreezing that produce backscatter variations could result in errors (Zhang, Wang, An, Liu, & Geng, 2022).Although we used a common strategy in considering the waveform parameters in the plane fit model to reduce this characteristic for Cryosat-2 data (Flament & Rémy, 2012;Sørensen et al., 2018;Zhang et al., 2020), the influence of the time-variable penetration depth might not be completely eliminated (Slater et al., 2019).
The impact of complex topography could result in slope-induced error.For Cryosat-2 data, its footprint size has an approximately 1.65 km footprint diameter.Such a large radar footprint cannot maintain track over regions with complex terrain and might cause the surface polynomial to approximate the local ice surface topography inaccurately (Brenner, Blndschadler, Thomas, & Zwally, 1983).We used large amounts of Cryosat-2 observations to fine-tune them for each grid cell.For laser data, it is not a concern because of the narrow divergence of the laser beam (Brenner, DiMarzio, & Zwally, 2007).The measurement error of ICESat and ICESat-2 can be as low as 0.1 m and 0.03 m, respectively (Shen, Ke, Yu, Cai, & Fan, 2021).Moreover, crossover analysis for laser data can reduce the elevation error caused by the slope (Felikson et al., 2017).
For ICESat/ICESat-2 data, the major errors in the laser retrieval are due to inaccurate modeling of the saturation and forward scattering effects on the return and the inaccuracies in the retrieved laser pointing angle.In addition, forward scattering caused by thin

Figure 1 .
Figure 1.Overview distribution of ICESat crossover points (a); number of Cryosat-2 data points in 5 km grid (b); and number of crossover points from ICESat-2 data in 5 km grid (c).

Figure 4 .
Figure 4. Volume change rates of the GrIS drainage basins divided by the elevation contour line of 2000 m from the three datasets.

Figure 5 .
Figure 5. (a) ATM L4 elevation changes in 2003-2009; (b) the difference in elevation change rates between ICESat and ATM L4; (c) ATM L4 elevation changes in 2010-2018; and (d) the difference in elevation change rates between CryoSat-2 and ATM L4.

Shuang
Liang received the Ph.D. degree in cartography and geographic information system from the Aerospace Information Research Institute (AIR), Chinese Academy of Sciences (CAS), Beijing, China.Her current and previous research interests include remote sensing of snow and ice.Huabing Huang received his Ph.D. degree from the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences.He is currently a professor in School of Geospatial Engineering and Science, Sun Yat-Sen University.His major research interests focus on dynamics land cover monitoring, vegetation height/biomass mapping with spaceborne Lidar.Xinwu Li is a professor of the Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing, China.He received a SAR Remote Sensing Ph.D. degree from Institute of Remote Sensing Application, CAS.His current research interests include Polarimetric and Interferometric SAR Remote Sensing.

Table 2 .
Elevation and volume change rates in each drainage system for three different periods.

Table 3 .
Schröder et al. (2019)lts of validation with ATM L4 observations.Schröder et al. (2019)over the Antarctic Ice Sheet.In the histogram, we can also see the distribution of differences present a normal distribution, and the data with a deviation between ±0.1 m/yr is the most distributed.Hence, our results are reliable.

Table 4 .
Comparison of total volume change rates for the GrIS derived by different studies.