A new Greenland digital elevation model derived from ICESat-2 during 2018-2019

. Greenland digital elevation models (DEMs) are indispensable to fieldwork, ice velocity calculation,s and mass 10 change estimations. Previous DEMs have provided reasonable estimations for the entire Greenland, but the time span of applied source data may lead to mass change estimation bias. To provide a DEM with a specific time-stamp, we applied approximately 5.8×10 8 ICESat-2 observations from November 2018 to November 2019 to generate a new DEM, including the ice sheet and glaciers in peripheral Greenland. A spatiotemporal model fit process was performed at 500 m, 1,2, and 5 km grid cells separately, and the final DEM is posted at the modal resolution of 500 m. 98% of the grids were obtained by 15 the model fit, and the rest DEM gaps were estimated via the ordinary Kriging interpolation method. Previous DEMs provided Greenland elevation information for different periods, but long temporal coverage introduced additional time uncertainty to scientific research. To provide a high-resolution DEM with a definite time, approximately 5.8×10 8 ICESat-2 observations from November 2018 to November 2019 were used to generate a new DEM for both the ice sheet and glaciers in peripheral Greenland. A spatiotemporal model fit process was first performed at 500 m resolution. To improve ICESat-2 20 data utilization, DEMs with 1 km and 2 km resolution across all of Greenland and an additional 5 km resolution in southernmost Greenland were used to fill the DEM gaps. Kriging interpolation was used to fill the remaining 2% of void grids that were insufficiently observed by ICESat-2 measurements. Compared with IceBridge mission data


Introduction
Greenland's digital elevation model (DEM) is particularly important for fieldwork planning and numerical modelling verificationThe digital elevation model (DEM) of Greenland is particularly important for fieldwork planning, numerical modelling verification, and ice movement tracking (Bamber et al., 2009;Bamber et al., 2013).The ice deformation rate and the underlying bedrock condition can be measured with the ice thickness data, which is useful to Combining the DEM with the ice thickness data, the ice deformation rate and the underlying bedrock condition can be measured, which can be used to determine subglacial hydrological pathways (Bamber et al., 2013).The surface elevation at different periods is also indispensable for studying elevation and mass changes to understand ice dynamics and estimate potential sea level changes (Sutterley et al., 2014;Smith et al., 2020).In addition, InSAR estimation of ice velocity requires high-accuracy and up-todate DEMs to distinguish phase differences caused by terrain and ice sheet movement (Riel et al., 2021).
The firstpreviously published Greenland DEM can dates back to the 1980s, providing elevations of peripheral Greenland generated through 3500 photographs from 1978 to 1987 with a resolution of 25 m, and 3500 photographs obtained from 1978 to 1987 were used to provide the elevations of peripheral Greenland with a resolution of 25 m (Korsgaard et al., 2016).
However, the low-visibility contrast between snow and ice surfaces may affect the radiometric and geometric quality of stereoscopic DEMs (Noh and Howat, 2015), which may introduce considerable uncertainty to the elevation.However, the data acquisition was limited by the low-visibility contrast between snow and ice surfaces (Noh and Howat, 2015), which introduced large time uncertainty into the DEM.Research regions were also restricted to the margin and outlets of Greenland, and there is a lack of understanding about the internal Greenland ice sheet.
The currently available DEMs of entire Greenland include those based on stereo-photogrammetry, altimeters, and interferometry.Most DEMs were derived from stereo-photogrammetry images, such as the Greenland Ice Mapping Project (GIMP) DEM (version 1) derived from ASTER, SPOT 5, and AVHRR photoclinometry (Howat et al., 2014), GIMP2 and ArcticDEM derived from GeoEye-1 and WorldView-1/2/3.ArcticDEM was the latest released DEM, with the highest resolution (2 m) among all free available Greenland DEMs.Optical image pairs may be influenced by weather, clouds, and the solar elevation angle (Korona et al., 2009); thus, the posted DEM is the combination of images of long timespan, which might limit its scientific applications into mass balance researches.In addition, owing to the wide coverage (86°N-86°S), high single-point accuracy (0.1-0.15 m), and small footprint size (70 m) (Zwally et al., 2002), ICESat has the ability to measure the elevation of entire Greenland.Hence, a bi-quadratic surface to fit ICESat footprints within each 1 km grid was adopted to obtain the ICESat DEM, but the largest search radius of 20 km in the low-latitude regions to some extent limited the ability to describe the small-scale elevation patterns at the Greenland margin (DiMarzio et al. 2007).A Ku-band synthetic aperture interferometric radar altimeter (SIRAL) carried by CryoSat-2 further increases the spatial coverage within 88°N-88°S.Although the footprint size (approximately 300 m) was larger than that of ICESat, the smaller cross-track distance (2.5 km) still ensures its ability to monitor the ice sheet (Wingham et al., 2002); thus, CryoSat-2 L1B data from 2011 to 2014 were applied to generate Cryosat-2 DEM through Kriging interpolation approach (Helm et al., 2014).Coarse across-track resolution is the major limitation to applying laser altimeters to generate a DEM with finer resolution (<1 km) in Greenland.TanDEM-X and TerraSAR-X were also used to generate the Greenland DEM using differential interferometry (Zink et al., 2014), but the X-band radar signal penetration depth into the dry snowpack may cause the elevation to be underestimated by several meters (Dehecq et al., 2016).These DEMs provide reasonable estimations for the entire Greenland, but the specific timestamps of the current DEMs were missing.
More DEM products became available to the public in the past 20 years with the improvement in observation methodologies.
ASTER and SPOT 5 measurements were used on the ice sheet periphery, and AVHRR photoclinometry was used in the inner and polar regions to derive the Greenland Ice Mapping Project (GIMP) DEM (version 1) (Howat et al., 2014).The GIMP1 DEM was then vertically calibrated by the Ice, Cloud, and Land Elevation Satellite (ICESat) and had a mean difference of -0.72 m with respect to the contemporaneous IceBridge data (Xing et al., 2020).ICESat had the advantages of wide coverage (86°N-86°S), high single-point accuracy (0.1-0.15 m), and small footprint size (70 m) compared with the former radar altimeter (e.g., Envisat, footprint size: 2-10 km, coverage: 81.5°N-81.5°S)(Zwally et al., 2002), enhancing the ability to measure the elevation of Greenland.Hence, the ICESat DEM adopted a bi-quadratic surface to fit all ICESat footprints within each 1 km grid, but the largest radius of 20 km in the low-latitude regions to some extent limited the ability to describe the small-scale elevation patterns at the Greenland margin (DiMarzio et al. 2007).
Benefitting from the ability to penetrate clouds, radar data can provide elevation information unaffected by weather conditions.CryoSat-2 carried a Ku-band synthetic aperture interferometric radar altimeter (SIRAL), further increasing the spatial coverage to 88°N-88°S.Although the footprint size (approximately 300 m) was larger than that of ICESat, the smaller cross-track distance (2.5 km) significantly improved its ability to monitor the ice sheet compared with that of ICESat (25 km) (Wingham et al., 2002).Thanks to the above advantages, CryoSat-2 L1B level data from 2011 to 2014 provided a reliable data source to provide elevation estimates (Helm et al., 2014).TanDEM-X and TerraSAR-X, high-resolution and allweather SAR data, were used to generate the Greenland DEM using differential interferometry (Zink et al., 2014), but the radar signal can penetrate into the snow, which causes the elevation to be underestimated.The GIMP2 elevation dataset was generated by the high-resolution panchromatic stereoscopic images of the GeoEye-1, WorldView-1, WorldView-2, and WorldView-3 satellites, and the time scale of the DEM was from 2009 to 2015.
ArcticDEM, the latest released DEM, which was generated by the same sources as GIMP2, has the highest resolution (2 m) among all free available polar DEMs.However, it is difficult to screen optical image pairs as DEM data sources because of the influence of weather, clouds, and the solar elevation angle (Korona et al., 2009).As a result, the mosaicked DEM is the combination of images from many times, and it is difficult to quantify the exact time of the DEM, limiting its scientific applications.
ICESat-2, a new generation of satellite-borne Lidar lidar altimeters, is intended as a successor to the ICESat mission to quantify the contribution of polar ice sheets to sea level rise and the impact of climate change (Markus et al., 2017).ICESat-2 has an orbital altitude of 500 km and an orbital inclination of 92°, accompanied by a revisit period of 91 days, which can provides centimetre-scale measurements of different surface types.The ICESat-2 beam footprint of approximately 17 m with a spatial interval of 0.7 m ensures accurate elevation measurements at a high orbital resolution by determining the local ice sheet slope (Neumann et al., 2019)The ICESat-2 beam footprint is approximately 17 m, with a spatial interval of 0.7 m (Neumann et al., 2019), which ensures accurate measurements of elevation at a high orbital resolution by determining the local ice sheet slope.A much finer observation can be obtained owing to its along-track distance of 0.7 m and cross-track distance of 3.3 km, which is a significant improvement compared with CryoSat-2's along-track distance of 0.3 km and crosstrack distance of 1.5 km, and ICESat's along-track distance of 170 mwhich is a great improvement over CryoSat-2's alongtrack distance of 1.5 km and cross-track distance of 3 km, and ICESat's along-track distance of 170 m.Not only the resolution but also the accuracy is has been improved.The accuracy in the flat ice sheet can reach 3 cm, and the accuracyit can still be less than 14 cm even for complex topography (Shen et al., 2021), which makes ICESat-2 a great data source to generate a DEM with high resolution and accuracy.
HenceHere, we present a novel Greenland DEM (ICESat-2 DEM) in May 2019 with a 500 m resolutionwe present a newly generated Greenland DEM with a resolution of 500 m using a spatiotemporal model fit, which is based on ICESat-2 measurements from November 2018 to November 2019.The overall accuracy of ICESat-2 DEM was evaluated by comparing it to the spatiotemporally matched IceBridge data.IceBridge data were used to evaluate the accuracy for all of Greenland and for different basins.The performance was also compared with other published DEMs under various terrain conditions to validate the reliability of the ICESat-2 DEM.

ICESat-2 ATL06 data
The ICESat-2 land ice height product ATL06 (Release 003) was used here for DEM generation.The product provides longitude, latitude, and surface heights based on the WGS84 ellipsoid and data acquisition time.The ATL06 product is developed from global geo-located photon data (ATL03) to estimate the land ice height (Smith et al., 2019).Compared with the original ATL03 product, land ice height is determined after instrument bias corrections (e.g., transmit pulse shape bias correction and first-photon bias correction) (Markus et al., 2017).The beam pair separation of the ATL06 product is set at 3.3 km across the track.The three pairs contain one strong beam and one weak beam, and the two beams within each pair are separated by 90 m distance.Brunt et al. (2019) compared the elevation of ICESat-2 ATL06 product and GPS data, and found that the accuracy differences of strong and weak beams are less than 2cm.Shen et al. (2021) compared ICESat-2 ATL06 product with IceBridge data under complex terrain, and the result indicated that the height difference between them is also trivial.Hence, we included weak beams to increase spatial coverage and data point utilization due to no systematic errors were found in strong and weak beams in ICESat-2 elevation measurements.The signal of the strong beams is stronger than that of the weak beams, so strong beams can provide more accurate measurements than weak beams.However, for strong and weak beams in the ATL06 product, both beams in one pair show similar performance, with a median difference of -0.08 cm and -0.13 cm for strong beam2 and weak beam1 (Shen et al., 2021).Although weak beams do not perform as well as strong beams, the accuracy is still higher than that of other altimeters.CryoSat-2, for example, has an accuracy of approximately 50 cm over the ice shelves and the interior of the ice sheet, with an error of more than 4 m on slopes greater than 0.9° (Wang et al., 2015).Hence, we did not exclude weak beams considering increased spatial coverage and data point utilization since no systemic error was found in ICESat-2 elevation measurements by different beams.However, only data marked as good quality (atl06_quality_summary=0) were used for DEM generation to improve the accuracy of the DEM.Over the entire Greenland ice sheet and ice capsoutlet glaciers, we used approximately 5.8×10 8 ICESat-2 elevation footprints to generate a new DEM, that is, the ICESat-2 DEM..

IceBridge data
To evaluate the accuracy of the ICESat-2 DEM, we used ATM surface elevation data from the IceBridge survey ATM surface elevation data from the IceBridge survey were used.ATM was intended to fill the gap between ICESat and ICESat-2, working at the same wavelength (532 nm) as ICESat-2.The absolute elevation accuracy of the ATM system can reach 0.1 m, and the position accuracy on the flat ice sheet is less than 1 m (Kurtz et al., 2013).The IceBridge ATM L2 Icessn Elevation, Slope, and Roughness Version 2 dataset was used to evaluate the DEMs.The final resolution of IceBridge was resampled to 25 m, and the estimated error was approximately 12 cm (Krabill et al., 2004).The root-mean-square error (RMSE) was taken as the roughness of each IceBridge data.The slope and aspect of IceBridge were described in Shen et al. (2021).The rootmean-square error (RMSE) was taken as the roughness of each IceBridge data.The slope and aspect were calculated as follows (Shen et al., 2021).
During 2009-2019, IceBridge provided millions of footprints over Greenland, which covering both the peripheral and inland areas of Greenland.The distribution of IceBridge data for May 2019, which was used to evaluate the accuracy of the new ICESat-2 DEM, is displayed in Figure 1.We also calculated the histogram of the elevation, surface slope, surface aspect, and roughness of IceBridge in May 2019.Overall, the elevations of the sampled regions ranged from 0 m to 3500 m, the surface slopes ranged from 0° to 10°, the surface aspects ranged from 0° to 360°, and the roughnesses ranged from 0 cm to 20 cm (Figure 2).These sampled areas had variable surface terrain conditions, which provided a reliable dataset to evaluate the performance of the ICESat-2 DEM.

Currently Other available other Greenland DEMs
We used other published DEMs to compare the performance of the generated ICESat-2 DEM, and the detailed information concerning these DEMs is provided in Table 1, and all DEMs have been referenced to the WGS84 ellipsoid.

GLAS/ICESat 1 km Laser Altimetry DEM
Greenland's DEM, derived from GLAS/ICESat laser altimetry data (from February 2003 to June 2005), provides the surface elevation for both Greenland ice sheets and caps, with less impact on slopes compared with radar altimetry data such as EnviSat and ERS 1/2.The spatial resolution is 1 km.The horizontal coordinates are based on polar stereographic coordinates system, and elevations are provided with respect to the WGS84 ellipsoid.

ArcticDEM
ArcticDEM is a high-resolution, high-quality digital surface model (DSM) of the Arctic in different spatial resolutions (2 m, 10 m, 32 m, 100 m, 500 m, and 1000 m) and its temporal coverage is from 2015 to 2018.ArcticDEM is a high-resolution, high-quality digital surface model (DSM) of the Arctic, with the highest resolution of 2 m, followed by resolutions of 10 m, 32 m, 100 m, 500 m, and 1000 m.The temporal coverage of the ArcticDEM project is mainly from 2015 to 2018.The mosaicked DEM files results are compiled from the best quality strip DEM files, and the filtered ICESat altimetry data are applied to improve the absolute accuracy.The estimated accuracy is approximately 85 cm at a resolution of 100 m (Xing et al., 2020).We used the elevation products of 500 m and 1000 m for comparisons.

TanDEM DEM
The TanDEM-X DEM (TanDEM) is a global DEM with a resolution of 90 m provided by the German Aerospace Centre (DLR).Data collection was completed in 2015, and global DEM production was completed in 2016 and published in 2018.
Different from previous datasets, it was generated by two X-band radar satellites (TanDEM-X and TerraSAR-X), which provides synchronous information to create a high-accuracy DEM aboutaccurate elevation information about the Earth's land surface.The absolute horizontal and vertical accuracy can reachis less than 10 m.The temporal coverage of the TanDEM data is mainly from 2011 to 2014.

CryoSat-2 DEM
CryoSat-2 L1B level data from January 2011 to January 2014 were used to provide the elevation of Greenland.Here, Helm et al. (2014) used waveform re-tracking to process LRM data and applied interferometry to process SARIn data, respectively.
Besides, he also leveraged slope correction to improve original product elevation accuracy.The biasCryoSat-2 L1B level data from January 2011 to January 2014 were used to provide the elevation of Greenland.Based on CryoSat-2 L1B level data from 2011 to 2014, waveform re-tracking and interferometry were performed on LRM and SARIn data, respectively, and slope correction was applied to the original product to improve the elevation estimates (Helm et al., 2014).The accuracy of the CryoSat-2 DEM was less than 1 m in flat regions and less than 4 m in rugged regions, showing similar performance to the other DEMs obtained by laser and radar altimeters.

DEM generation
We followed the method of Slater et al. (2018) to compute the elevation of Greenland, which is an iterative least-squares fit model to all the elevation measurements in each grid as a quadratic surface.This model is described by Equation (1).To compute the elevation of Greenland, we followed the method of Slater et al. (2018), which is a spatiotemporal model to fit the surface height in each grid.The model used is described as follows: where hi is the elevations derived from ICESat-2 measurement points in one grid, h represents the modelled elevation, and a0, a1, a2, a3, and a4 are surface elevation fluctuations.dh/dt is the elevation change rate in the 13 months, t is the month difference between May 2019 and ICESat-2 acquization time.tmid is the time of the mid timestamp (May 2019), and (x,y) are the coordinates in the polar stereographic projectionwhere hi represents the modelled elevation, dh/dt is the elevation change rate in the 13 months, t is the ith month starting from November 2018, tmid is the time of the mid timestamp (May 2019), and (x,y) are the coordinates in the polar stereographic projection.The uncertainty is the 95% confidence level for elevation.
We found that there are more voids in low-altitude areas due to the low density of ICESat-2 footprints during the procedure of Greenland DEM generation.Therefore, it is necessary to select a suitable spatial resolution to accomplish DEM with fewer gaps.DEM accumulated within 250 m gridsThe selection criterion for DEM resolution is to depict the detailed elevation pattern, so finer resolution should be adopted.A much finer resolution can reveal more detailed elevation patterns; this will lead to more calculated gaps, making the proportion of interpolated grids larger and leading to larger biases and uncertainties.The proportions of the calculated grids are latitude-dependent, and high-latitude areas account for higher coverage of the calculated grids, owing to their dense orbital distribution.However, it is difficult to increase the proportion of calculated grids in low-latitude areas by using high spatial resolution.
Hence, we tested resolutions of 250 m, 500 m, 1 km, 2 km, and 5 km to achieve the optimal resolution.We found that the 250 m grid covered only covered 15.38% of the Greenland area and only 30% even at high latitudes, ; hence,so we discarded this resolution for further processing.In contrast, a 500 m resolution increases overall coverage to 33% and nearly 70% at high latitudescoverage at high latitudes to nearly 70% (Figure 3).With a 1 km and 2 km resolution, the proportion of calculated portions grids exceeds 90% in the basins regions that are north of 75°N (basins 1, 2, and 8).The resolution of 2 km covers 82% of all Greenland.However, a 2 km resolution cannot obtain optimal coverage in low-elevation areas, while a 5 km resolution can further increase the coverage of calculated grids to 98%, especially in the southern basins (basins 4, 5, and 6).
To reduce the effect of any poor fit, quality control was constrained in terms of data availability, quality, and rationality.We set the minimum number of grid points to 10 and the minimum timestamp to 2 months, which could ensure that enough measurements were contained in a grid cell to generate a reliable DEM.We assumed that the maximum elevation change is 10 m/yr and its uncertainty is impossible to exceed 0.4 m/yr (Slater et al., 2018).Furthermore, we assumed that DEM uncertainty is less than 10 m and the maximum RMSE in each grid is 10 m.After this filter procedure, the elevation range is -500 m-3600 m and this result is feasible since it is within the elevation range of published Greenland DEM products.
After the aforementioned process, we have acquired Greenland DEM in four resolutions (500 m, 1 km, 2 km, and 5 km).
However, these four types of DEM all include void areas thus we need to incorporate them to obtain final Greenland DEM results with the minimal gaps.Firstly, we used Greenland DEM with 500 m resolution as our first DEM source.Afterwards, Greenland DEMs with 1 km, 2 km, and 5 km resolution were resampled to 500 m by applying a bilinear method to fill the gaps in this DEM and the finer resolution as our first option.Unavoidably, there are still some voids in the final Greenland DEM, but this has a minor impact on DEM accuracy.In this study, we described the unvoided area (98%) in the final Greenland DEM as 'calculated grids' and termed the rest (2%) as 'interpolated grids'.For the rest, an ordinary kriging approach was used to interpolate.The ICESat-2 DEM was posted at the modal resolution of 500 m after gap filling and interpolation.A median filter of 2.5 km×2.5 km was applied to the posted ICESat-2 DEM to minimize the influence of different resolutions.
Here, we applied different methods to estimate ICESat-2 DEM uncertainties in calculated grids and interpolated grids.The elevation uncertainty of the calculated grids was calculated by Equation ( 2) based on MATLAB R2018a.For interpolated grid uncertainty estimation, we just used kriging variance error calculated by ArcGIS 10.6.There is a 95.5 percent probability that the actual elevation at the grid is the predicted raster value ± two times the square root of the variance error of the corresponding cell by assuming the kriging errors are normally distributed.Hence, the two times the square root of the value in the variance error was taken as the elevation uncertainty in the interpolated grids (Equation ( 3)).
elevation uncertainty interpolated grids = 2 × √variance error where b  is the elevation, SE (b) is the standard error of the elevation, and t (1-0.025,n-p) is the 95% percentile of tdistribution with n-p degrees of freedom, n is the number of ICESat-2 measurements in one grid, p is the number of regression coefficients (i.e., 7), and variance error is kriging variance error.
The slope was calculated by the method of Horn et al. (1994), thus the slope uncertainty was calculated based on the law of propagation: where   is the elevation uncertainties of the adjacent grids of the central grid.
Ultimately, we separated Greenland into 8.5×10 6 grids to compute the elevation in a 500 m×500 m grid.To build a fineresolution DEM and decrease the interpolation coverage as much as possible, the 1 km and 2 km results were used to fill the observed gaps in the DEM by resampling their results to 500 m by bilinear interpolation.The same process was repeated at a 5 km resolution in southernmost Greenland.The results of each resolution were constrained in terms of data availability, data quality, and rationality.We set the minimum number of grid points to 10 and the minimum timestamp to 2 months, which could ensure that enough measurements were contained in a grid cell to make a good elevation estimation.In addition, we introduced thresholds to remove outliers, which are RMSE≥10 m, the uncertainty of elevation change ≥10 m, the rate of elevation change≥10 m/yr and the rate of elevation change uncertainty ≥0.4 m/yr (Slater et al., 2018).The original kriging method was adopted to fill the 2% data holes that still existed.To reduce the influence of different resolutions, a median filter of 2.5 km×2.5 km was applied to the final ICESat-2 DEM.

DEM accuracy evaluation
One ICESat-2 DEM grid cell usually has several IceBridge measurement points.In each grid cell, the ICESat-2 DEM elevation values were subtracted from the median of all IceBridge elevations within it, and this difference was seen as the final bias.First, we grouped IceBridge data according to different DEM grids and calculated the median elevation, slope, aspect, and roughness within each grid, taken as the true values in this grid.The height differences were obtained by subtracting the DEM elevation at the location of IceBridge from the median IceBridge elevation.Subsequently, we used the median difference (MED), the mean difference (MD), the median absolute difference (MAD), the standard deviation (STD), the RMSE, the 90% confidence interval LE90 (Gonzalez-Moradas and Viveen, 2020; Höhle and Höhle, 2009), and the correlation (R) to evaluate each DEM.The calculations are as follows: where ℎ  is the elevation difference in each DEM grid and n is the number of overlapping IceBridge footprints.
We additionally used the elevation intervals of 0 m to 500 m, 500 m to 1000 m, 1000 m to 1500 m, 1500 m to 2000 m, and ≥2000 m to study the relationship between the elevation difference and elevation.For the surface slope, we divided the slope into 5 intervals of 0° to 0.25°, 0.25° to 0.5°, 0.5° to 1°, 1° to 2°, and ≥2° to detect the relationship between the elevation difference and slope.Similarly, the same step was repeated for roughness intervals of 0 cm to 5 cm, 5 cm to 10 cm, 10 cm to 15 cm, 15 cm to 20 cm, and ≥20 cm.We identified the aspect as north, east, south, and west to investigate the relationship between the elevation difference and terrain aspect.
measurements.Hence, the other grids were interpolated from the neighbouring elevations by using the ordinary kriging method.
The ICESat-2 DEM shows the same pattern as the other published Greenland DEMs, .The highest elevation appears in the ice sheets and shows a downward trend to the margins (Figure 4 (a)), and large topographic fluctuations occur on the outlet glaciers around the periphery of Greenland.Furthermore, the monthly elevation change rate was also obtained from spatiotemporal model fit, thus the DEM for each month from November 2018 to November 2019 can be derived theoretically.with the highest elevation appearing in the interior ice sheets and showing a downward trend to the margins.
Lower elevations occur on the surrounding individual glaciers around the periphery of Greenland, which also exhibits large topographic fluctuations.In addition to the DEM, the model fitting method can also be used to obtain the monthly elevation change rate; hence, the DEM for each month from November 2018 to November 2019 can be derived theoretically.
The DEM uncertainty and slope uncertainty show obvious latitude-dependent patterns.Larger values tend to be found at low latitudes, and this pattern may be related to the number of ICESat-2 measurement points in each grid cell.The uncertainty also presents an increasing trend from the interior to the margins, which is approximately less than 0.5 m in the inner ice sheet and higher uncertainty of 2-5 m can be observed for the periphery of Greenland (Figure 4 (b)).The generated slope uncertainty is large at the edges, which is also concurrent with slopes exceeding 1° (Figure 4(c) and 4(d)).The accuracy of satellite laser altimeters is affected by surface roughness, slope, and other environmental factors (Brunt et al., 2017).A flatter surface provides a more uniform reflection than a steeper surface of the measurement footprint, and the more accurate height measurements of the original ICESat-2 footprints can be obtained in the low-slope regions, hence the higher accuracy of the ICESat-2 DEMWe calculated the mean elevation of the calculated grids at four resolutions according to different main basins (Figure 5).It can be clearly seen that with increasing resolution, the calculated elevations were mostly higher, while the calculated elevation were decreased only in basin 1 under 2 km resolution and in the GLA region under 1 km and 2 km resolution.The regions with the largest bias were concentrated in the low-latitude regions of Greenland, with a maximum deviation of more than 40 m.The elevation difference did not increase with varying resolutions for the GLA region.Because the topography here was more complex, it was possible that the elevation could be either overestimated or underestimated on different glaciers, and this uncertainty alleviated the elevation differences.

Evaluation of ICESat-2 DEM by comparing with IceBridge data
The ICESat-2 DEM compares favourably to those of IceBridge data (Figure 5 and Table 2).The bias between Greenland DEM and IceBridge data in calculated grids is smaller than that of interpolated grids, which indicates that the elevations derived from model fit tend to be more accurate than those estimated from interpolation.Poorer performance in the interpolated grids is reasonable due to the low spatial correlation in the regions with large surface fluctuations like the Greenland south margins.The application of three resolutions may add additional effects, i.e., different grid cell resolutions tend to present different elevation estimates, thus DEMs of different resolutions (namely, 500 m, 1 km, 2 km) were evaluated (Table 3).The results show that the DEM with 2-km resolution exhibits worse performances than those of 500m and 1km resolutions, but the bias of 2 km-resolution DEM is still smaller than the error in the interpolated grid.It is reasonable to apply three resolutions for posting ICESat-2 DEM of 500 m resolution considering the differences between the calculated and interpolated grids.
We also compared the accuracies according to the 10 basins covered by data from May 2019 (Table 4), the accuracy of the ICESat-2 DEM shows an apparent spatial trend that better accuracy is observed in the north than in the south basins, and the pattern may be related to the small proportion of calculated grids in the southern basin and the application of DEMs with 2 km and 5 km resolution.We calculated the mean elevation of main basins at four resolutions to further assess the effects of different spatial resolutions when generating ICESat-2 DEM (Figure 6).It can be seen that the calculated elevations were generally higher with increasing resolution, and the regions with the largest bias were concentrated in the low-latitude basins of Greenland.The small elevation difference for the GLA region was possibly caused by the elevation being overestimated or underestimated on different glaciers due to the complex topography, and this uncertainty alleviated the elevation differences.
There are still some differences between the ICESat-2 DEM and IceBridge data, mainly due to the inconsistent coverages of the two datasets.It should be noted that the IceBridge data used for evaluation were distributed only at latitudes below 75°N, where the posted DEM was mostly derived from DEMs at coarse resolutions.The ICESat-2 DEM should have higher accuracy in the regions located beyond the north of 75°N.Hence, these biases are acceptable because the evaluated value represents the upper bound of the ICESat-2 DEM bias, and the deviation should be smaller when considering ICESat-2 DEM as a whole.

ICESat-2 DEM uncertainty
The uncertainty in the DEM shows an obvious spatial pattern.The uncertainty presents an increasing trend from the interior to the margins and from the north to the south.In the inner ice sheet, the uncertainty is approximately less than 0.5 m, but for the periphery of Greenland, higher uncertainty can be observed (Figure 4

Comparison with other available DEMs
We first compared the spatial coverage of different DEMs and found that the coverage percentage of all DEMs exceeded 99% (Table 1).Of these, the CryoSat-2 DEM has the lowest coverage (99.33%), and most peripheral glaciers show a large number of voids.TanDEM has the second lowest coverage (99.93%), and although the release product has filled the gaps, there are still small void grids over peripheral Greenland due to radar shadow, echo delays and phase unwrapping errors.The ICESat DEM has the largest coverage (99.99%), which is due to its 5.5-20 km search radii for different latitudes when 带格式的: 标题 1 calculating the elevations of each grid.The ICESat-2 DEM coverage is comparable to that of the 500 m and 1 km ArcticDEM (99.98%).
The elevation differences between the new ICESat-2 DEM and the other five published DEMs show that DEMs usually perform better for low-slope regions (Figure 7).The ICESat-2 DEM is generally close to that of the 500 m ArcticDEM except in complex terrains, which can prove the great reliability of the ICESat-2 DEM.In particular, significant positive values can be seen in the elevation difference between the ICESat-2 DEM and TanDEM on the Greenland ice sheet, which is assumed to be caused by the X-band penetration into the snowpack.All the difference maps show significant negative values in the Jakobshavn Isbrae glacier, where experienced the greatest loss of the Greenland ice sheet (Smith et al., 2020), this phenomenon may reflect the real elevation changes during different DEM acquisition times.
Figure 7 depicts the elevation differences between the new ICESat-2 DEM and the other five published available DEMs.
Larger elevation differences can be found in peripheral Greenland than those in all other DEMs, but the elevation differences in the interior ice sheets are much smaller.Significant positive patterns can be seen in the elevation difference between the ICESat-2 DEM and TanDEM on the Greenland ice sheet due to X-band penetration.All the difference maps show significant negative values in the Jakobshavn Isbrae glacier, which is the area experiencing the greatest loss of the Greenland ice sheet (Smith et al., 2020).This reflects the real elevation changes in Greenland at different DEM acquisition times.The ICESat-2 DEM is generally close to that of the 500 m ArcticDEM except in complex terrains.As satellite optical images have a much finer original spatial resolution (2 m) than ICESat-2, larger elevation differences in high-slope regions are reasonable.However, these differences are much smaller than those of other DEMs, which proves the greater reliability of the ICESat-2 DEM.
In this study, we used only IceBridge data that overlapped with the corresponding DEM period to evaluate all DEMs' vertical accuracies for the entire Greenland and areas with little elevation changes (-0.05~0.05m/yr) (Smith et al, 2020).Results show that the ICESat-2 DEM showed significant improvements in accuracy compared with other altimeter-derived DEMs, and is also comparable to DEMs derived from stereo-photogrammetry and interferometry.Compared with ICESat DEM derived from the 6.9×10 6 footprints and CryoSat-2 DEM derived from the 7.5×10 6 footprints, approximately 80 times as many data were used to generate the DEM, thus the finer-resolution of 500 m and higher accuracies can be obtained (Table 5 and 6).The ICESat-2 DEM has the best performance with regard to the parameters except for MED of CryoSat-2 DEM.
Contrary to expectations, there is no elevation underestimation in the CryoSat-2 DEM, possibly because slope and topographic corrections have been performed.TanDEM has less bias in the coastal region, this is possibly caused by that ICESat data have been used to calibrate the raw DEMs there (Wessel et al., 2016).Although model-based or empirical models can, to some extent, correct the penetration bias (Abdullahi et al., 2019), such correction in the ice sheet is generally restricted to the regional scope (Wessel et al., 2021), thus significant elevation underestimations cannot be corrected and still exist.The accuracy of the ICESat-2 DEM is higher than the 1 km ArcticDEM, and is comparable to the 500 m ArcticDEM.
ICESat data were used to calibrate the ArcticDEM in both the horizontal and vertical directions to increase the accuracy of ArcticDEM, but the actual changes in the ice surface introduce additional uncertainties because ICESat data predate the accuracy (Candela, 2017).
The median differences in surface slope and roughness for these DEMs illustrate that all their elevation biases become larger with increasing slope and roughness with the exception of TanDEM (Figure 8(b) and 8(c)).500 m ArcticDEM should maintain high accuracy as stereo-photogrammetry can provide more consistent elevation values at the regional scale than altimeter, the similar performance of ICESat-2 DEM can verify its stability under all terrain conditions.Larger elevation differences in high-slope regions are reasonable as satellite images have a much finer original spatial resolution (2 m) than ICESat-2.Furthermore, DEMs generated by stereo pairs have obvious directivity in terms of surface aspect (Figure 8(d)).
The accuracy on the north slope is significantly lower, mainly due to the poor illumination condition of the images in the north direction.
We evaluated the vertical accuracies of all DEMs using IceBridge data as described in section 3.2.To ensure comparability between datasets, we used only IceBridge data that overlapped with the corresponding DEM period for evaluation.
Compared with those of DEMs generated by altimeters, the accuracy of the DEM is greatly improved (Table 4).The ICESat-2 DEM generally performs the best.For all of Greenland, the Cryosat-2 DEM has the smallest MED of 0.03 m, and the ICESat-2 DEM ranks next.The ICESat-2 DEM has the best performance with regard to the rest of the parameters.The correlation can reach 0.9999 for all of Greenland.The ICESat-2 DEM shows the best performance in areas with elevations greater than 2000 m.In areas with elevations less than 2000 m, the ICESat-2 DEM still shows the best performance, except that the smallest MED appears in the Cryosat-2 DEM.Contrary to expectations, there is no elevation underestimation in the CryoSat-2 DEM, possibly because slope and topographic corrections were performed in the DEM.
Compared with DEMs generated by image pairs, the ICESat-2 DEM also shows good performance (Table 4).Affected by the penetration of SAR signals into dry snow in high-elevation areas, TanDEM has the worst performance in the regions with elevations above 2000 m, corresponding to an MED of -3.76 m.The MED (-2.32 m) is lower in the regions with elevations below 2000 m, indicating that TanDEM is less affected by penetration in the low-elevation regions.The snow penetration of the X-band is mainly concentrated in dry snow areas, and the penetration into wet snow in low-elevation areas is limited and can even be neglected.In addition, ICESat data alone were used to calibrate the raw DEMs in the coastal regions (Wessel et al., 2016), and the effect of penetration was alleviated to a certain extent.All parameters indicate that the accuracy of the ICESat-2 DEM is higher than that of the 1 km ArcticDEM.The accuracy is highly comparable to that of the 500 m ArcticDEM.The MED is better for all of Greenland, and the rest of the parameters are similar.However, the ICESat-2 DEM has a higher elevation correlation with IceBridge, especially in regions with elevations greater than 2000 m.
The accuracy of all DEMs is affected by varying topography.All DEMs tend to increase in accuracy with increasing elevation, with the exception of TanDEM (Figure 8 .The accuracy on the north slope is significantly lower.This is mainly due to the poor illumination condition of the images in the north direction, which affects the accuracy of the generated DEM.The accuracy of satellite laser altimeters is affected by surface roughness, slope, and other environmental factors (Brunt et al., 2017).The measurement footprint is sensitive to changes in surface terrain, and a flatter surface provides a more uniform reflection than a steeper surface.A more accurate height measurement of the original ICESat-2 footprints can be obtained in the low-slope regions.
In summary, the ICESat-2 DEM is superior to the previous satellite altimeter-derived DEMs in both spatial resolution and elevation accuracy.Compared with the 6.9×10 6 footprints used in the ICESat DEM and the 7.5×10 6 footprints used in the CryoSat-2 DEM, approximately 80 times as many data were used to generate the DEM; thus, the ICESat-2 DEM with higher resolution of 500 m can be obtained.Due to the characteristics of ICESat-2 itself, the denser orbit density and higher measurement precision, higher elevation accuracy can be achieved.Smaller differences can be found when compared to the satellite image-derived DEMs.When comparing the DEMs generated by SAR interferometry, the ICESat-2 DEM has an increased accuracy because it is not affected by the penetration depth into snow.The penetration depth depends on the radar wavelength and the snow and ice properties, leading to a significant underestimation of elevation (Guerreiro et al., 2016).
Although model-based or empirical models can correct the penetration bias to some extent (Abdullahi et al., 2019), such correction is generally restricted to the regional scope (Wessel et al., 2021), and there are still significant elevation underestimations.The ArcticDEM has high precision because ICESat data were used to calibrate the ArcticDEM in both the horizontal and vertical directions to increase the accuracy.However, the ICESat data predate the ArcticDEM by almost 10 years, and actual changes in the ice surface introduce additional uncertainties.Finally, there may also be systematic errors between the ArcticDEM's different sensors (Candela, 2017).
The problem remains with the uncertainty the ignorance of the elevation changes during DEM data acquisition.The generation method of the ICESat-2 DEM is similar to that of the ICESat DEM, but the advantage of the former is that we took the data acquisition time of each point into account and used the smaller search radius in the interpolation, minimizing the time uncertainty as much as possible

Conclusions
A new digital elevation model of Greenland was provided based on the ICESat-2 observations acquired from November 2018 to November 2019.The DEM was posted at a modal resolution of 500 m.98% of the grids were directly derived from a model-fit method, and the further 2% were interpolated by kriging method.The application of different resolutions can reduce the number of interpolated grids and less bias in the elevation estimation.Compared with spatiotemporally matched elevation measurements from the IceBridge dataA model-fit method was applied within the grid cells within different spatial resolutions to estimate the surface elevations with a modal resolution of 500 m.Different resolutions ensured that more ICESat-2 data could be utilized to reduce the number of interpolated grids, which contributed to less bias in the elevation estimation.Overall, we estimated the uncertainty with a median difference of -0.48 m for the entireall of Greenland, which represents the upper bound of the ICESat-2 DEM bias.The accuracy of the ICESat-2 DEM shows an apparent spatial trend, and better accuracy can be observed in the northern than in the southern basins owing to the denser coverage of ICESat-2 tracks in the high-latitude regions.
In addition, the accuracy of the ICESat-2 DEM shows an apparent spatial trend, and better accuracy can be observed in the north than in the southern basins, owing to the denser coverage of ICESat-2 tracks in the high-latitude regions.
Compared with other published Greenland DEMs, i.e., the ICESat DEM, CryoSat-2 DEM, 1 km ArcticDEM, 500 m ArcticDEM, and TanDEM, the ICESat-2 DEM maintains great accuracy stability under various topographic conditions.The ICESat-2 DEM is superior to the previous satellite altimeter-derived DEMs in both spatial resolution and elevation accuracy.
(b)).The glaciers in northern Greenland exhibit uncertainties of 2-5 m, while the uncertainty increases to 10 m in southern Greenland.The generated slope and slope uncertainty are in good agreement with the elevation uncertainty, and the uncertainty is large at the edges, which is also concurrent with slopes exceeding 1° (Figure4(c) and (d)).The error for the independent glaciers is due to the large fluctuation in the surface slope, and the ICESat-2 DEM with a resolution of 500 m has difficulty capturing detailed ICESat-2 DEM and IceBridge show general agreement (Figure6and Table2).Overall, the entire area of Greenland exhibits MED, MD, and MAD values of -0.48 m, -1.90 m, and 2.73 m, respectively, which compare favourably to those of IceBridge data.The MEDs do not differ in the calculated and interpolated grids (both -0.48 m).However, the performance decreases in the interpolated regions, which can be seen from the remaining parameters, such as MD, MAD and STD.The correlations between ICESat-2 DEM and IceBridge are high for all of Greenland, the calculated grids, and the interpolated grids, with the highest correlation occurring in all of Greenland and the calculated grids (both 0.9999).Kriging interpolation is used to fill the remaining 2% of the DEM gaps, which are concentrated in the southernmost part of Greenland.Although ICESat-2 has a higher coverage percentage than previous altimeters, interpolation errors caused by the data gaps cannot be avoided.Poorer performance in the interpolated grids is reasonable due to the low spatial correlation in the regions with large surface fluctuations.Although this accuracy is worse than that of the calculated grids, it still agrees well with that of the IceBridge data.We also compared the accuracies according to the different basins.The data from May 2019 covered only 10 basins in Greenland (Table3).Regionally, the MED of each basin ranges from -0.95 m to 0.20 m, with the largest deviation appearing in basin 5.The MAD, STD, RMSE, and LE90 show the same trend, with the largest bias in the GLA region.The highest MAD is located in basin 5, excluding the GLA region.For the calculated grids, all basins display a difference of less than 5 m, but the difference in the GLA region can be up to 15 m.For the interpolated grids, the poorest performance still appears in basin 5 and the GLA region.Overall, the accuracy of the ICESat-2 DEM shows an apparent spatial trend.Better accuracy is observed in the north than in the south, mainly because the proportion of calculated grids in the southern basin is small, and the results of coarse resolution are mostly adopted.Additionally, due to higher temperature sensitivity, true elevation changes during the acquisition times of IceBridge and ICESat-2 lead to another part of the bias.The ICESat-2 DEM has very high accuracy, but there are still some differences between the ICESat-2 DEM and IceBridge data, which are mainly due to the inconsistent resolutions of the two datasets.It should be noted that the IceBridge data used for evaluation were distributed only at latitudes below 75°N.The DEM was mostly derived by the results with a 1 km resolution or even with 2 km and 5 km resolutions.The ICESat-2 DEM should have higher accuracy in the regions with 500 m resolution (located north of 75°N).Hence, these biases are acceptable because the evaluated value represents the upper bound of the ICESat-2 DEM bias, and the deviation should be smaller when considering Greenland as a whole.
(a)).The ICESat-2 DEM is clearly superior in regions with elevations less than 500 m compared with DEMs generated by other altimeters.It is also superior to the 1 km ArcticDEM at elevations below 2000 m, and the accuracy is close to that of the 500 m ArcticDEM.Similarly, elevation differences between each DEM and IceBridge were calculated with respect to the surface slope (Figure8(b)).Similar to the elevation results, only TanDEM shows a trend of improvement in accuracy as the slope increases and has the highest accuracy in the high-slope (also the low-elevation) regions.The CryoSat-2 DEM, 1 km ArcticDEM and ICESat DEM all tend to show higher errors as the slope increases, while both the ICESat-2 DEM and 500 m ArcticDEM maintain high accuracy and stability under all slope conditions.A trend of decreasing performance is found with increasing roughness for all DEMs except TanDEM (Figure8(c)).The ICESat-2 DEM and 500 m ArcticDEM have the best performance, with deviations of less than -0.5 m for all roughnesses.In terms of aspect, DEMs generated by stereo pairs have obvious directivity(Figure 8(d)) Smaller elevation differences between ICESat-2 DEM and DEMs derived from stereo-photogrammetry and interferometry can imply the reliability of the ICESat-2 DEM.Compared to DEMs derived from satellite altimeters, the ICESat-2 DEM has the finest spatial resolution for a Greenland DEM and hence has more accurate elevation measurements.Obvious elevation underestimation can be seen in the TanDEM due to penetration into snow.Compared to the ArcticDEM, the ICESat-2 DEM has smaller differences in the Greenland ice sheet, which implies the reliability of the ICESat-2 DEM.Although the uncertainties in the ICESat-2 DEM are affected by the ICESat-2 measurements themselves, and the DEM in the low-latitude regions were derived from results of coarse spatial resolution, and the coarse spatial resolution in the low-latitude regions, the specific time-stamped ICESat-2 DEM can benefit studies of elevation change and mass balance in Greenland.theadvantage of accurate time can benefit studies of elevation change and mass balance in Greenland.More ICESat-2 data can be used to generate DEMs with higher resolution as more ICESat-2 observations become availableWith more ICESat-2 observations becoming available, more ICESat-2 data can be used to generate DEMs with higher resolution, especially in the southernmost glaciers.

Figure 1 :Figure 2 :
Figure 1: IceBridge data acquired in May 2019, which were used to evaluate the generated ICESat-2 DEM, covering regions in Greenland with various terrain conditions: (a) elevation, (b) slope, (c) aspect, and (d) roughness.Labels in the picture are the main glaciers in Greenland.640

Figure 4 :
Figure 4: (a) Elevation of the Greenland DEM calculated from 13 months of ICESat-2 footprints acquired between November 2018 and November 2019.(b) The elevation uncertainties.(c) Slope derived from the elevation map (a) and (d) Slope uncertainty 650

Figure 5 :
Figure 5: Elevation differences of 9 main regions under different resolutions, which are calculated by subtracting the 500 m DEM 655

Figure 6 :
Figure 6: Elevation difference calculated as IceBridge data subtracted from the new ICESat-2 DEM.(a) Calculated grids and (b) interpolated grids.IceBridge data were acquired in May 2019.

Figure 7 :
Figure 7: Elevation differences calculated between the new ICESat-2 DEM and the other five published available DEMs.(a) ICESat DEM, (b) CryoSat-2 DEM, (c) 500 m ArcticDEM, (d) 1 km ArcticDEM, and (e) TanDEM.For each picture, the previously published DEM was resampled to 500 m, and the difference was calculated as the resampled DEM subtracted from the new 665

Figure 8 :
Figure 8: Elevation differences between different DEMs and IceBridge data under different terrain conditions.(a) Elevation, (b) slope, (c) roughness, and (d) aspect.The solid black lines near the box centres denote median values of elevation differences, the upper and lower boundaries of each box denote upper and lower quartiles (Q1 and Q3), the length means the interquartile range 670