Edinburgh Research Explorer Assessment of terrain elevation estimates from ICESat-2 and GEDI spaceborne LiDAR missions across different land cover and forest types

,


Introduction
With the launch of two new National Aeronautics and Space Administration (NASA) spaceborne Light Detection and Ranging (LiDAR) missions in late 2018 -the Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) and the Global Ecosystem Dynamics Investigation (GEDI) -a new suite of three-dimensional (3D) data of Earth's land surface has been produced. These missions are designed to provide accurate vertical information about land surfaces, such as vertical vegetation structure, and terrain elevation. Information on vertical vegetation structure is required to assess the status and temporal dynamics of forests and provides baseline data to estimate terrestrial carbon stocks and emissions (Baccini et al., 2012;Saatchi et al., 2011). Information on terrain elevation is crucial for applications such as hydrological modeling (Sanders 2007;Wilson et al., 2007) and modeling of soil erosion (De Vente et al., 2013). Furthermore, accurate ground detection is key for estimation of vertical vegetation structure parameters from LiDAR (e.g., relative vegetation heights, and vegetation density at different height ranges), since these structure metrics are calculated as height above the ground. Terrain elevation data below the forest canopy can also aid in estimating vegetation height from global Digital Elevation Models (DEMs) (e.g., Kellndorfer et al., 2004;Simard et al., 2006). For this, however, the penetration depth of the signal must be considered (Schlund et al., 2019).
LiDARs are active systems that send pulses of light and measure their travel time. Some systems also measure the returned intensity. Knowing precisely the travel time of the signal, the distance to scattering objects can be calculated (Lim et al., 2003). LiDARs for terrestrial applications usually operate in the near-infrared (NIR) range of the electromagnetic spectrum between 900 and 1545 nm. Due to its ability to resolve small gaps within tree crowns, LiDAR is able to detect the ground under dense canopies as well as derive vegetation structure and understory information. A LiDAR system illuminates an area (the instantaneous laser footprint) by sending a laser pulse, which can generate one or many returns depending on the Earth's surface characteristics (bare-ground or vegetated). The laser footprint can vary in size from centimeters for an airborne sensor to tens of meters for a spaceborne instrument (e.g., 25 m for the GEDI instrument).
ICESat-2 was launched in September 2018 and carries an innovative instrument, the Advanced Topographic Laser Altimeter System (ATLAS). In contrast to its predecessor ICESat, the ICESat-2 mission follows a different acquisition concept and uses different technology. Unlike the waveform LiDAR sensor used by ICESat GLAS (Geoscience Laser Altimeter System), ATLAS is a photon-counting LiDAR (PCL) that uses less energy per pulse, allowing high repetition rates (Markus et al., 2017). High repetition rates improve sampling density along track, collecting laser footprints every 70 cm (Markus et al., 2017). In contrast to the dual-wavelength green (532 nm) and infrared (1064 nm) ICESat GLAS instrument, ICESat-2 ATLAS operates 532 nm (green) lasers with a repetition rate of 10 kHz (Markus et al., 2017). The choice of wavelength was set by the low dead-time and efficiency of the single photon detectors. The illuminated footprint size of ATLAS is around 14-17 m (Markus et al., 2017;Popescu et al., 2018) at the beginning of the mission and will increase to around 20 m at the end of the 3-year planned mission (Neuenschwander and Pitts 2019). ICESat-2 possesses three laser beam pairs, which are separated by ca. 3.3 km (Neuenschwander and Pitts 2019), allowing data collection along six parallel tracks and thus improving spatial coverage. Each beam pair comprises a strong and weak beam with an energy ratio of 4:1 and a spatial separation of 90 m, which helps to determine cross-track slope (Markus et al., 2017;Popescu et al., 2018). The main application of ICESat-2 is monitoring of the cryosphere, i.e., mapping of changes in glacier and ice-sheet elevation and corresponding contributions to sea-level change and estimation of sea ice thickness. Additionally, ICESat-2 collects LiDAR samples over land areas and provides height measurements over global forests. The operation of a laser with a wavelength of 532 nm is optimal for ice and snow monitoring. However, this wavelength is not ideal for vegetation monitoring due to a reduced reflectance from leaves compared with near-infrared, higher background solar noise (Swatantran et al., 2016) and lower atmospheric transmittance.
In December 2018, the GEDI sensor was successfully launched and installed on the International Space Station (ISS). Because the GEDI lasers operate from the ISS, the data collection coverage is limited to regions between 51.6 • N and 51.6 • S, thus covering temperate and tropical forests but omitting most boreal forests. GEDI was specifically designed for vegetation monitoring and operates three 1064 nm lasers (Dubayah et al., 2020a). Two are full-power lasers while the third laser is split into two coverage beams, generating a total of four beams. Coverage beams use approx. 1/3 the energy of the full-power beams (Dubayah et al., 2020a), resulting in reduced capability to receive ground returns under dense canopies . To improve spatial coverage, the four beams are dithered on alternative shots, producing eight parallel ground tracks (Dubayah et al., 2020a). The tracks are separated by around 600 m across-track and 60 m along-track. GEDI is a full-waveform LiDAR with a footprint size of around 25 m. After the nominal two-year mission by the end of April 2021, GEDI had generated more than 10 billion cloud-free observations (Dubayah et al., 2020a). Fortunately, the GEDI mission is extended until September 2023, allowing collection of many more observations.
In this study, we assessed the accuracy and precision of terrain elevation estimates from these two new spaceborne LiDAR missions, since terrain elevation errors are the primary drivers of other product errors (e.g., vegetation height, plant area index, aboveground biomass) Hancock et al., 2019). For this, we examined the "ATLAS/ICESat-2 L3A Land and Vegetation Height", version 5 (Neuenschwander and Pitts 2019; and the "GEDI Level 2A Footprint Elevation and Height Metrics", version 2 (Dubayah et al., 2020b). Uncertainties in spaceborne LiDAR-based terrain estimates can be related to atmospheric conditions (haze), ranging precision of the instrument, horizontal geolocation accuracy, terrain relief, seasonal effects (e.g., leaf-on/-off canopy conditions, inundation), and algorithms applied to determine ground return. The accuracy of the ground detection is primarily a function of the amplitude of the ground return above noise, which is driven by canopy cover.
We compared the performance of terrain elevation estimates from the spaceborne LiDAR data in different biomes, using high-resolution airborne measurements as reference data. We performed our analysis over six regions spanning a range of different land cover classes, vegetation types, canopy heights, tree densities, and levels of anthropogenic disturbance. Instead of assessing terrain elevation estimates in different geographical sites, we analyzed the accuracy of spaceborne LiDAR terrain estimates across land cover and forest types. For this, we merged the data from different sites and assigned them to a certain land cover and forest type. Finally, we investigated the impact of different environmental conditions (terrain slope, canopy height, and canopy cover), and sensor parameters (e.g., beam sensitivity (for GEDI) and number of terrain photons (for ICESat-2)) on the accuracy of spaceborne LiDAR products. The objective of the study is, however, not to compare performance between the spaceborne sensors (due to different footprint size), rather to provide statistics of current available versions under different conditions (land cover, canopy cover, terrain slope).

Study areas
We investigated the accuracy of both sensors in four different land cover classes and six forest types located in Brazil, Germany, South Africa and the USA (Fig. 1). The first two sites are located in temperate broad-leaved and needle-leaved forests in central Germany (50-51 • N). The third study site includes mixed land cover types in Sonoma County, California, USA (38-39 • N). The fourth and fifth sites are located in the Brazilian Amazon (2.5-10 • S): one site is spatially distributed over six areas of pasturelands and secondary forests, and one is located within intact tropical rainforest in the mid-Juruá region. The sixth study site is situated in a southern African savanna ecosystem in Kruger National Park, South Africa (23-25 • S).

Land cover and forest types
Since the accuracy of terrain elevation estimates relates to land cover and forest types, we merged data from different sites and assessed their accuracy at two levels. At the highest level we analyzed the accuracy of spaceborne LiDAR data for four land cover classes, which cover the major part of our study sites: bare soil, herbaceous vegetation, savannas and forests. We masked out urban areas in Sonoma (USA), Hainich and Roda (Germany), as well as tidal wetlands in Sonoma and open water for Amazon sites. Furthermore, we did not assess the accuracy of spaceborne LiDAR terrain estimates for agricultural areas, since this class includes vineyards in Sonoma and orchards in Hainich and Roda sites, and represents a mixture of herbaceous and woody vegetation.
Land cover information was derived from high-resolution land cover maps for Sonoma (SonomaVegMap 2019, 2020), Hainich and Roda (TLBG 2020a). These land cover maps are based on aerial imagery, partly on machine learning techniques (for the Sonoma site) and fine-scale manual editing (SonomaVegMap 2019; TLBG 2020a). In our analysis, bare soil areas are natural non-water, barren or sparsely vegetated areas. The herbaceous vegetation class consists of non-woody vegetation in Sonoma, Hainich and Roda, as well as areas with canopy height (derived from airborne laser scanning (ALS) data) lower than 2 m in the Amazon sites. All ICESat-2 and GEDI samples collected over Kruger National Park were assigned to the savanna land cover class. Forested areas include pixels classified as forest in Sonoma (Sonoma-VegMap 2019, 2020), Hainich and Roda (TLBG 2020a) and areas with an ALS canopy height greater than 5 m in the Amazon sites.
At the lower level we subdivided the forest land cover class into forest types occurring in our study areas: temperate broad-leaved, temperate needle-leaved, temperate mixed, tropical secondary, tropical floodplain, and tropical upland forests. A detailed description of forest types is given in the next subsection.
Temperate needle-leaved forests are represented by Scots pine (Pinus sylvestris), Norway spruce (Picea abies) and rarely occurring European larch (Larix decidua) from Hainich and Roda, and coast redwood (Sequoia sempervirens) and Douglas fir (Pseudotsuga menziesii) from the Sonoma site.
Mixed forest stands with the abovementioned tree species are assigned to the temperate mixed forests class.
Based on ALS data, the temperate needle-leaved forest class possesses the tallest trees (mean canopy height 21.1 m) and densest tree canopies (mean canopy cover 78.9%) among the temperate forest classes (Fig. 2). Canopy cover (CC) is defined as the percentage of all ALS point clouds from heights greater than 5 m. In contrast, temperate broad-leaved forests feature smaller mean canopy height (14.5 m) and lower canopy cover (64.3%) compared to other temperate forest classes (Fig. 2).

Tropical forests
Classification into tropical intact and tropical secondary forests was accomplished using geographical locations of the ALS tracks and "The Intact Forest Landscapes" (IFL) map for the year 2016 (Potapov et al., 2017), and confirmed by observation of historical optical satellite imagery. Twelve ALS tracks in the Juruá study region (Amazonas State) and within the binary IFL map were assigned to the class intact tropical forests. Six airborne LiDAR transects, which were acquired across the states of Acre, Pará and Roraima, are located outside the binary IFL map, and with ALS canopy height higher than 5 m were classified as tropical secondary forests. These six ALS tracks were collected over a mixture of Terrain slope for four land cover classes (a) and forest types (b); canopy height (c) and canopy cover (d) for six forest types. The parameters are calculated from ALS point cloud data or (for Kruger National Park) Digital Mapping Camera (DMC) stereoscopic data at a pixel size of 25 m. Canopy Height (RH95) is defined as the 95th percentile height of point cloud data. Canopy cover (CC) is defined as the percentage of all returns from heights greater than 5 m. TemBF -temperate broad-leaved forests; TemMF -temperate mixed forests; TemNF -temperate needle-leaved forests; TroBFF -tropical broad-leaved floodplain forests; TroBUF -tropical broad-leaved upland forests; TroSFtropical broad-leaved secondary forests. We further separated tropical intact forests of the mid-Juruá region into upland ("terra-firme") and floodplain forests using a Synthetic Aperture Radar (SAR)-based Amazon-wide wetland map (Hess et al., 2015); upland/floodplain boundaries were refined using recent high-resolution Landsat and PALSAR-2 imagery. Uplands and floodplains of the mid-Juruá region support largely undisturbed, biodiverse forests (Albernaz et al., 2012;Cardoso et al., 2017;Wittmann et al., 2006). Floodplain forests include those of the Juruá River (a high-sediment, "whitewater" tributary of the Amazon River), as well as smaller floodplains of local tributaries with lower sediment loads. About half of the Juruá floodplain forest occupies flats and depressions inundated from 3 to 10 months per year with seasonal water depths up to several meters, and half is situated on higher ground (levees, scroll bar ridges) with flooding depths of less than 2 m. Hawes et al. (2012) reported stand densities of 633 and 639 stems per hectare for floodplain and terra-firme forest respectively, for the mid-Juruá region. Based on 12 ALS tracks, mid-Juruá terra-firme forests are taller and have greater canopy cover than floodplain forests (Fig. 2).

Savanna
The Kruger National Park (KNP) represents the savanna study site. KNP is the largest national conservation area in South Africa with an extent of almost 19,500 km 2 . Vegetation is characterized by a mixture of growth forms, primarily grasses, shrubs, and trees. In large parts of the area trees are deciduous and shed their leaves during the dry season (May-September), while both grass and woody vegetation are green during the wet season (October-April). The dominant tree species are Mopane (Colophospermum mopane) in the dry northern part (mean annual precipitation (MAP) below 500 mm/year), Red bushwillow (Combretum apiculatum) and Knob Thorn (Acacia Nigrescens) in the more humid southern part of the park with MAP above 500 mm/year (Eckhardt et al., 2000;Venter et al., 2003). According to a woody cover map based on C-band SAR data (Urban et al., 2020), mean woody cover of the KNP is 23%. Thus, the KNP site can be considered as an open savanna with high woody cover (up to 80%) along the riparian areas (Venter et al., 2003). Since the KNP is a conservation area, the changes in woody cover are induced primarily by large herbivores such as elephants (Asner and Levick 2012) or by fires. Based on terrain height estimates derived from aerial imagery, the elevation in the park varies between 105 and Fig. 4. Absolute terrain height difference (GEDI terrain height (default parameter "elevation lowestmode" varies among six algorithms depending on plant functional type) minus reference terrain height) plotted as a function of GEDI degrade flag (a), absolute deviation between GEDI terrain and TanDEM-X DEM heights (b), GEDI time of acquisition (c), GEDI number of modes (d), GEDI beams (e), GEDI beam sensitivity (f). 832 m with a mean elevation of 333 m and a mean slope of 3.4 • (Heckel et al., 2021). Mean vegetation height from a high-resolution canopy height model is around 5 m with the tallest (95th percentile) canopies reaching 12 m.

Reference data
As reference data we used high resolution airborne data. For Hainich, Roda, Sonoma, and the Brazilian Amazon the data were collected with LiDAR sensors providing detailed three-dimensional information. In the Kruger National Park reference information on ground topography was derived from aerial stereoscopy using an optical sensor (color-infrared channels). Terrain slope, canopy height and canopy cover estimates grouped by land cover and forest types based on the reference data are shown in Fig. 2. A detailed description of each reference data set is given in the next subsections.

Hainich and Roda ALS
Small-footprint discrete-return airborne LiDAR data in Hainich were collected wall-to-wall by the Thuringian land surveying office (TLBG) between December 2016 and March 2017. The mean point density is around 20 returns per m 2 . ALS data at the Roda site were collected for the entire area by the TLBG during two flight campaigns: about 80% of the data were acquired in February 2014, and the remaining 20% were recorded in November 2014 and March 2015. The mean point density of this ALS data is ca. 10 returns per m 2 . The original data are freely available at (TLBG 2020b) in three formats: gridded digital terrain model (DTM) (bare ground) and digital surface model (DSM) (including object heights, e.g., trees and buildings) at 1 m pixel spacing as well as point clouds in.laz format. All data possess a horizontal coordinate reference system (CRS) ETRS89, UTM Zone 32 N and GRS80 ellipsoid, and a vertical reference system DHHN2016 with heights above a German Combined QuasiGeoid 2016 (GCG2016). The reported horizontal accuracy of the ALS data is 15 cm and vertical accuracy is 9 cm. We converted the horizontal datum to Latitude/Longitude in the WGS84 horizontal reference system (EPSG: 4326). We further transformed vertical datum to the EGM2008 geoid using the GCG2016 geoid (i.e., first, from GCG2016 orthometric heights to ellipsoidal heights; then from ellipsoidal heights to EGM2008 orthometric heights).
As a reference for bare-earth elevation we used gridded DTM at 1 m pixel spacing, and the canopy height model (CHM) was derived by subtracting DTM from DSM. ALS data were collected during the winter season under leaf-off conditions.

Sonoma ALS
Wall-to-wall airborne LiDAR data over Sonoma County, California were collected between September and November 2013. Reported mean pulse density of the data is 11 pulses per m 2 (WSI 2016). The Sonoma ALS data are publicly available in the form of point clouds (.las format), gridded terrain height, canopy height, and canopy cover layers at 1 m pixel spacing (SonomaVegMap 2020). Based on 9348 ground check points, maximum deviation of the ALS elevation is 16 cm, with an RMSE of 3 cm (WSI 2016).
The horizontal datum of the data is NAD 1983 StatePlane California II FIPS 0402. The vertical datum of ALS data is North American Vertical Datum of 1988 (NAVD88) based on GEOID12A. The units of original data are in US survey feet. We converted the horizontal datum to Latitude/Longitude in the WGS84 horizontal reference system (EPSG: 4326), and the vertical datum to the EGM2008 geoid. Further, US survey feet were recalculated in meters.

Amazon pasture and Juruá forest ALS
Small-footprint discrete-return airborne LiDAR data for the two Amazonian sites were acquired by Brazil's National Institute for Space Research (INPE) from April to July in 2016 as part of its Amazon Biomass Estimation (EBA) project (Dalagnol et al., 2021;Tejada et al., 2019). In total, 18 of the 610 EBA transects were analyzed covering an area of ca. 6,900 ha. The mean point density in the ALS tracks is around 4 returns per m 2 . The pulse footprint was set to be below 30 cm, with horizontal accuracy within 100 cm and vertical accuracy within 50 cm. The original data are available as point clouds in.las format. All data use the SIRGAS2000 horizontal coordinate reference system and GRS80 ellipsoid. The vertical reference system is the Imbituba tide gauge station. The mean difference between the Imbituba tide gauge station and the EGM2008 geoid is around − 0.5 m (Gruber et al., 2012). Our results for non-forested areas (i.e., for the samples with canopy height smaller than 5 m derived from ALS data) indicate a mean difference between the Imbituba tide gauge station and the EGM2008 geoid to be around − 0.79 m, 0.15 m and 0.11 m for GEDI footprints, ICESat-2 100 m segments and ICESat-2 20 m segments, respectively (Fig. S7).
We pre-processed point clouds with LAStools software (version 200,304) to generate bare-ground elevation and vegetation metrics. On the point clouds we first applied the noise removal function "lasnoise" implemented in LAStools. We then classified points into ground and non-ground returns using the LAStools function "lasground". Classified ground returns were then triangulated and rasterized using the LAStools function "las2dem". The final ALS-based bare-ground elevation has a pixel spacing of 1 m. To generate vegetation metrics (canopy height (RH95) and canopy cover), we used the LAStools function "lascanopy".

KNP Digital Mapping Camera
To create very high-resolution elevation models for the entire Kruger National Park (KNP), Digital Mapping Camera (DMC) aerial imagery from the National Geo-spatial Information (NGI) programme of South Africa's Department of Rural Development and Land Reform (DRDLR) was utilized. Data sets were acquired from an altitude of approximately 5,500 to 6,000 m at a nominal pixel resolution of 25 cm in the time period between September and October 2018 (dry season). Aerial stereoscopic algorithms to derive height information (surface and terrain height) were carried out using the Enterprise software package from CATALYST (formerly known as PCI Geomatics), setting the vertical datum as the EGM2008 geoid. Major processing steps included meta data preparation, bundle adjustments and tie point collection, semiglobal matching (SGM)-based height extraction as well as several filtering steps as described in (Heckel et al., 2021). Validation of the DMC-based DTM with GNSS point samples yielded an LE90 (Equation (4), Section 2.5) of 1.02 m (Heckel et al., 2021). A canopy height model was produced by subtracting the DTM from the DSM. The retrieved DTM and CHM have a pixel spacing of 5 m.

Spaceborne LiDAR data
A short description of the spaceborne LiDAR data as well as product levels and versions used in this study is given in the following subsections. The total number of samples per land cover and forest type of the examined spaceborne LiDAR data is given in boxplots (Figs. 5, 7 and 8).

GEDI data
All available GEDI L2A version 2 data (Dubayah et al., 2020b) acquired over study regions between April 2019 and August 2021 were used for the accuracy assessment of terrain elevation estimates. GEDI L2A samples provide latitude, longitude, ground elevation, relative height metrics above the ground (as a proxy for canopy height) and surface energy metrics at footprint level (25 m), as well as TanDEM-X DEM height. The terrain elevation estimates of the L2A products are based on six different algorithms. These algorithms define thresholds for signal start and end, as well as smoothing width for noise and signal (for more details see Table 5 in Hofton et al. (2019)). In contrast to the GEDI version 1 data, where default terrain elevation estimates are based on algorithm 1, in version 2 the default terrain elevation estimates vary among these six algorithms depending on plant functional type (Beck et al., 2021). For instance, in areas with dense undergrowth, algorithms with a lower signal end threshold might be preferred.
In addition to the GEDI L2A data, we extracted EGM2008 geoid heights for each footprint from the L1B product (Dubayah et al., 2020c), in order to convert GEDI terrain elevation, which is referenced to the WGS84 ellipsoid, to orthometric height (relative to the EGM2008 geoid). Furthermore, estimated canopy cover and beam sensitivity parameters for each GEDI sample were retrieved from the L2B product (Dubayah et al., 2020d). Beam sensitivity is defined as the maximum canopy cover through which the GEDI can detect the ground with 90% probability . Before the accuracy assessment, we applied the following filters on the L2A data: 1) we selected GEDI samples with a quality flag of 1 (valid waveform) in all six algorithms; 2) we removed GEDI samples for which the absolute difference between estimated terrain elevation and TanDEM-X DEM was greater than 50 m for all six algorithms; 3) we removed samples for which estimated canopy cover provided in the GEDI product is greater than the corresponding beam sensitivity; 4) we selected non-degraded GEDI samples (i.e., GEDI degrade flag parameter equal to 0). Degraded conditions in GEDI L2A data represent, e.g., degraded attitude and degraded trajectory (Beck et al., 2021). Applying the GEDI degrade flag reduces the number of GEDI observations by 23%, from 1.78 million to 1.37 million (Fig. 4a). As our GEDI terrain elevation estimate we used the parameter "elev_lowestmode", which represents the elevation of the center of the lowest mode relative to the WGS84 ellipsoid. We used only those GEDI samples that were located completely within the reference data.  5. Absolute terrain height difference (GEDI terrain height (default parameter "elevation lowestmode") minus reference terrain height) plotted as a function of land cover (a), forest types (b), ALS terrain slope (c), ALS terrain slope for non-forested areas (ALS CHM <5 m) (d), ALS canopy height (e), ALS canopy height for flat terrain (with ALS slope <5 • ) (f), ALS canopy cover (g), ALS canopy cover for flat terrain (with ALS slope <5 • ) (h). Note: there are a lower number of samples on impact of ALS canopy cover, since this information is not available for the savanna site. The boxes show the interquartile range (25th and 75th percentile), the bold bar the median and the vertical bars are 1.5 times the interquartile range; outliers are excluded for improved visualization. Number for each slope group gives corresponding number of samples. TemBFtemperate broad-leaved forests; TemMFtemperate mixed forests; TemNF -temperate needle-leaved forests; TroBFF -tropical broad-leaved floodplain forests; TroBUFtropical broad-leaved upland forests; TroSFtropical broad-leaved secondary forests.

ICESat-2 data
For the accuracy assessment of ICESat-2 data we used the ATL08 product (ATLAS/ICESat-2 L3A Land and Vegetation Height) version 5  collected between October 2018 and September 2021. The first step of the ATL08 algorithm is the classification of raw photons into noise and signal photons. The ATL08 product is based on geolocated photons from the ATL03 product (Neumann et al., 2020), which contains a signal flag indicating a likelihood that a photon is classified as a signal. As reported in Neuenschwander and Pitts (2019), however, this default algorithm for the ATL03 product works over smooth surfaces such as ice sheets, but sometimes classifies photons from the top of canopy as noise. To improve the detection of signal photons for the vegetated surfaces, an additional signal finding algorithm DRAGANN (Differential, Regressive, and Gaussian Adaptive Nearest Neighbor) is applied (Neuenschwander and Pitts 2019). A more detailed description of the ATL08 product and corresponding photon classification from the ATLAS sensor can be found in (Neuenschwander and Pitts 2019;. It is expected to receive up to four photons per single strong beam footprint over vegetated surfaces, which is a much lower number than over highly reflective ice and snow surfaces (ca. 10 signal photons) (Neuenschwander and Pitts 2019). In order to have enough ground and canopy photons, and thus, provide more robust estimates, ATL08 canopy and terrain elevations are estimated for transects of 100 m length, i. e., 140 sequential footprints with a diameter of ca. 14 m are combined. Furthermore, ATL08 product version 5 provides estimates on ground elevation and canopy height for 20 m segments (i.e., ~28 sequential footprints are combined). We examined ATL08 product at both spatial resolutions (100 m and 20 m, hereafter ATL08_100 and ATL08_20).
For our ICESat-2 terrain elevation estimate, we used the parameter "h_te_best", which provides the best fit terrain elevation at the mid-point location of each 100 m and 20 m segments. We applied the following filters on the ATL08_100: we removed samples for which the absolute difference between estimated terrain elevation ("h_te_best") and the MERIT-DEM (Multi-Error-Removed Improved-Terrain DEM (Yamazaki et al., 2017)), provided in ATL08 data for each ICESat-2 transect, was greater than 30 m. The same absolute height difference threshold was applied by Neuenschwander et al., (2022). Further, we removed those ATL08_100 samples with estimated terrain uncertainty higher than 20 m. Determinations of both thresholds are based on comparison with reference data (Fig. 6). Finally, we selected only those ATL08_100 samples with maximum elevation height of 2,000 m and with maximum canopy height ("h_canopy") of 100 m. For the ATL08_20 product we applied a threshold of 50 m for an absolute height difference between estimated terrain elevation ("h_te_best") and the CopernicusDEM (Airbus 2020), since for 20 m segments neither MERIT DEM nor terrain uncertainty estimates are provided. This threshold is higher than that used for MERIT-DEM, since CopernicusDEM provides the elevation of the scattering height of the X-band signal. We used only those ICESat-2 samples that were located completely within the reference data.

Data extraction and calculation of statistics
To compare spaceborne-based terrain elevation estimates with reference data we performed the following steps. First, since the elevations in the reference data are orthometric heights referenced to the EGM2008 geoid (Sections 2.3.1-2.3.4), a conversion from the WGS84 ellipsoid to the EGM2008 geoid as vertical reference for GEDI and ICESat-2 data was applied. The EGM2008 geoid heights (Pavlis et al., 2012) are provided in the GEDI L1B product as auxiliary information. For ICESat-2 segments we extracted the EGM2008 geoid heights from the National Geospatial-Intelligence Agency 2.5' grid (Pavlis et al., 2012), since the geoid heights are not provided for the 20 m ATL08 product. Additionally, we masked out urban areas in Hainich, Roda and Sonoma, as well as tidal wetlands in Sonoma and open water for Amazon sites. Fig. 3 illustrates an example of an ALS-based CHM collected over tropical upland forests (a), vertical distribution of ALS point clouds over an ICESat-2 transect (b) and over a GEDI sample (d) together with corresponding GEDI waveform (c), as well as the main processing and analysis steps (e).
For terrain elevation comparison, we calculated a median value from reference data over GEDI and ICESat-2 footprints. For statistical assessment of terrain height estimates, the following metrics between spaceborne and reference data were calculated: Fig. 7. Absolute terrain height difference (ICESat-2 terrain height (parameter "h_te_best") at 100 m (ATL08_100) minus reference terrain height) plotted as a function of land cover (a), forest types (b), ALS terrain slope (c), ALS terrain slope for non-forested areas (ALS CHM <5 m) (d), ALS canopy height (e), ALS canopy height for flat terrain (ALS slope <5 • ) (f), ALS canopy cover (g), ALS canopy cover for flat terrain (ALS slope <5 • ) (h). Note: there are a lower number of samples for estimating impact of ALS canopy cover, since this information is not available for the savanna site. There are no ICESat-2 samples available for a canopy cover range of 90-95%.

Mean error
Median absolute deviation MAD = median(|Δh i − median(Δh i )|) Linear error 90% LE90 = quantile90(|Δh i |) where h i is a height estimate for a footprint from spaceborne data and h ref is the estimated median height within the corresponding footprint from reference data. These metrics have previously been used for accuracy assessments of global DEMs (Baade and Schmullius 2016; Höhle and Höhle 2009;Wessel et al., 2018). We did not report coefficient of determination (R 2 ), since this metric is dependent on local topographic variability and not the measurement accuracy. For instance, for a site with topography varying by 1,000 m, an uncertainty of tens of meters can still have an R 2 higher than 0.9, while a flat site will have a very low R 2 , even with sub meter accuracy.

Results
A summary of goodness-of-fit statistics of the GEDI and ICESat-2 terrain elevation estimates for different land cover and forest types is shown in Tables 1 and 2. Corresponding scatterplots are shown in the Section "Supplementary materials". In the following subsections results on impact of environmental conditions and sensor parameters on accuracy and precision of GEDI and ICESat-2 terrain estimates are described.

Data screening using GEDI auxiliary metrics
Auxiliary information available for every GEDI sample can help to Fig. 8. Absolute terrain height difference (ICESat-2 terrain height (parameter "h_te_best") at 20 m (ATL08_20) minus reference terrain height) plotted as a function of land cover (a), forest type (b), ALS terrain slope (c), ALS terrain slope for non-forested areas (ALS CHM <5 m) (d), ALS canopy height (e), ALS canopy height for flat terrain (ALS slope <5 • ) (f), ALS canopy cover (g), ALS canopy cover for flat terrain (ALS slope <5 • ) (h). Note: there are a lower number of samples for estimating the impact of ALS canopy cover, since this information is not available for the savanna site. select confident terrain estimates, i.e., those with a high probability of being true ground elevation. Outliers in GEDI terrain estimates can be caused by degraded pointing and positioning information, and reflection from clouds. Furthermore, vertical complexity, acquisition times and beam mode (power or coverage) can affect accuracy of terrain height estimates (e.g., Duncanson et al., 2020).
Applying the GEDI parameter "degrade_flag" helps to increase precision of elevation estimates by filtering out more than 20% of GEDI observations (Fig. 4a). Considering all land cover classes, the LE90 (as an indicator for precision) (Equation (4)) for non-degraded samples was 4.1 m, versus 12.9 m for degraded samples and 5.6 m for all samples (i. e., degraded and non-degraded). In the next steps, we investigated nondegraded GEDI samples only. Elevation height from global DEMs can further improve accuracy and precision of GEDI estimates (Fig. 4b).
Here, however, using a low absolute difference threshold (e.g., 10-20 m or 20-30 m) can lead to removing samples that were collected over forested areas, bearing in mind that TanDEM-X DEM provides scattering height of the X-band signal. Global DEMs and this absolute difference threshold can be used to screen GEDI samples acquired over clouds (e.g., by using a threshold of 50 m).
Although the proportion of degraded samples was lower for nighttime acquisitions (18%) than for day-time acquisitions (28%), we did not observe a higher accuracy and precision for night-time GEDI samples than for those collected during the day (Fig. 4c). Mean error was 0.59 m for night-time and 0.54 m for day-time samples. One reason for this might be that degraded day-time GEDI data had already been filtered out in previous steps. About 20% more GEDI samples were collected during the night than during the day. The GEDI parameter "number of modes" provides an estimated number of peaks in a waveform and reflects the vertical complexity of a GEDI sample. As expected, an increasing "number of modes" leads to a higher deviation between GEDI and reference data (Fig. 4d). GEDI beams 0 to 3 are so-called "coverage" beams with a reduced capability to receive ground returns under dense canopies . In our study regions, lower precision for beams 0 to 3 was found only for degraded samples; for non-degraded samples, beams 0 to 3 show similar correlation to reference data as the "power" beams (beams 5 to 11), though with a smaller number of samples (Fig. 4e). Since the GEDI "sensitivity" parameter provides an estimation of maximum canopy cover through which GEDI can detect the ground with 90% probability, this parameter should support selection of samples over dense forest, i.e., with higher sensitivity parameters denser canopy covers can be penetrated. However, in our study regions, a higher GEDI sensitivity parameter led to greater GEDI residual errors, which conflicts with the definition of this parameter (Fig. 4f).

Impact of environmental conditions on GEDI terrain elevation estimates
Accuracy of GEDI terrain elevation estimates varies between land cover and forest types ( Fig. 5a and b) and is a function of environmental conditions such as terrain slope, vegetation structure and canopy cover (Fig. 5c-h). Mean error of GEDI terrain elevation estimates across our land cover classes is below 1 m (Table 1). Thus, we can conclude that the GEDI default algorithm for estimation of terrain elevation works in most cases. The most precise GEDI terrain estimates (in terms of LE90) are found for land cover classes bare soil and savanna (Table 1, Fig. 5a).

Table 1
Goodness-of-fit statistics between spaceborne-based terrain elevation estimates and reference data sorted by land cover; "n" provides corresponding number of samples for each category.  Table 2 Goodness-of-fit statistics between spaceborne-based terrain elevation estimates and reference data sorted by forest types; "n" provides corresponding number of samples for each category. Most parts of the savanna landscape are characterized by scattered trees, so that in most cases GEDI has a strong ground return. Furthermore, because Kruger National Park is rather a flat area with a mean slope of 3.4 • (Fig. 2a), the impact of terrain slope is low for this site. Finally, due to a large number of samples, the statistics for this site are less affected by outliers. Terrain elevation estimates for herbaceous vegetation have a similar strong correlation with reference data as for bare soil and savanna. However, the precision metric (LE90) here is around twice that over bare soil and savanna (Table 1), meaning that the precision itself is lower. The lowest precision across land cover classes was found for forested areas (Table 1, Fig. 5a). This land cover class features the highest mean slope (Fig. 2a) and possesses a more complex vertical structure compared to other land cover types. Across the forest types, GEDI terrain estimates over temperate broadleaved, mixed and needle-leaved forests possess similar statistics with reference data (Table 1, Fig. 5b) with slightly lower RMSE and LE90 in temperate broad-leaved forests (Table 1). One reason for this can be caused by lower canopy cover and slightly flatter terrain compared to other temperate forest types ( Fig. 2b and d). Moreover, partly leaf-off conditions (winter 2019-2020 and winter 2020-2021) in temperate broad-leaved forests might improve detection of ground return.
In contrast to temperate forests and savanna, the agreement in terrain estimates from GEDI and ALS in the tropical broad-leaved forest was found to be much poorer, with an overestimation of terrain elevation for GEDI's default algorithm (Table 1, Fig. 5b). Dense canopies with multiple vegetation layers cause a weak return from the ground (e.g., Fig. 3c), which complicates an accurate classification of the terrain. In this environment a lower signal end threshold (e.g., GEDI algorithm 5) can improve correlation with reference data (Fig. 9b and e).
Considering all land cover and forest types together, increasing slopes lead to an increasing negative deviation between GEDI and reference terrain estimates (Fig. 5c). Furthermore, this increased deviation is also observed over non-forested areas (i.e., areas with a canopy height of lower or equal than 5 m) (Fig. 5d). The reason for this, is partly owing to horizontal geolocation error of the current GEDI version 2 data.
The mean geolocation accuracy of this data is 10 m (1 σ offset) (Beck et al., 2021). A horizontal offset of the footprint center of 10 m would translate to a 10 m*tan(slope) vertical offset. For instance, for a 30 • slope a mean vertical error of about 6 m would be expected. Growing residual errors of GEDI terrain estimates are also noted for canopy heights greater than 10 m for all terrain types (Fig. 5e). However, no correlation between the residual error and canopy height over flat terrains (slope <5 • ) was found (Fig. 5f). That is, canopy height itself did not impact the accuracy of ground detection. Similarly, denser canopies lead to a higher deviation of GEDI terrain estimates for all terrain types (Fig. 5g). Over flat terrains (slope <5 • ), only very dense canopies (>90%) cause high residual errors of GEDI terrain estimates (Fig. 5h).

Data screening using ICESat-2 auxiliary metrics
ICESat-2 samples contain auxiliary information that can help to screen outliers. We did not note differences between ICESat-2 beam modes (strong vs. weak) for terrain estimation (Fig. 6a). However, as there are more than twice as many samples collected by strong beam than those collected by weak beam, samples collected by weak beam could have been filtered out in the first pre-processing steps (Section Fig. 9. GEDI terrain elevation estimates based on default algorithm (upper panel) and alternative algorithm 5 (lower panel) for three tropical broad-leaved forest types: floodplain (left), upland (center), and secondary (right) forests. 2.4.2). Although elevation accuracy for ICESat-2 is expected to be lower during the day than at night owing to solar noise, we found slightly higher errors for night samples (Fig. 6b). A possible explanation is the filtering of lower-confidence day-time samples from the ATL03 product (which is a basis for the ATL08 product). As for GEDI's terrain residual error, an absolute deviation between ICESat-2 terrain and reference height estimates can be used to remove the outliers. Furthermore, ICESat-2 parameters "number of classified terrain photons" and "estimated terrain uncertainty" can further aid in refining ICESat-2 terrain estimates (Fig. 6c-f). Obviously, the more classified terrain photons there are, the more robust terrain estimation should be. Terrain uncertainty parameters comprise systematic uncertainties (e.g., pointing, geolocation) and errors from photon classification (e.g., noise, terrain, canopy) normalized by the number of terrain photons . Depending on how many ICESat-2 samples remain and what accuracy is needed, one can set a specific threshold using this auxiliary information.

Impact of environmental conditions on ICESat-2 terrain elevation estimates
In contrast to GEDI estimates, the highest mean error (ME) of ICESat-2 terrain estimates were found for bare soil and herbaceous vegetation with negative offsets of − 1.27 m and − 0.97 m, respectively (Table 1).
Similar to the GEDI assessment, ICESat-2 bare-ground elevation estimates showed the highest precision (in term of LE90) in the savanna site (Table 1, Fig. 7a). A likely explanation for the better ICESat-2 performance over savanna vs. bare ground is the approximately 500-fold greater sample sizes for savanna; because of the small numbers of samples for bare ground, a relatively small number of outliers (partly located on slopes) had a significant impact on accuracy and precision. The statistical metrics between ICESat-2 and ALS data across temperate forest types were similar (Table 1, Fig. 7b). The correlation between ICESat-2 and reference data was the lowest for the tropical broad-leaved upland forest compared to other forest types (Fig. 7b), but this is based on a very small number of samples (e.g., Fig. S4e).
Similar to the GEDI terrain height residual errors, we observed an increasing negative error of ICESat-2 terrain elevation estimates with increasing terrain slope for all samples (Fig. 7c) as well as for the samples over non-forested areas (Fig. 7d). Tall trees on flat terrain did not impact the accuracy of the ICESat-2 terrain elevation (Fig. 7f), similar to the negligible effect on GEDI terrain estimates. Over our study sites only dense canopies (>95%) impact ICESat-2 terrain elevation estimates (Fig. 7h).
In general, accuracy assessment of the ATL08_20 product showed similar trends as for the ATL08_100 product, though with a higher number of samples (Tables 1 and 2 ,Figs. 7 and 8). Accuracy of the ATL08_20 product for land cover classes and forest types followed a similar pattern the accuracy of the ALT08_100 product, i.e., it depends on vertical complexity of illuminated areas (e.g., multiple vegetation layers).
We observed an increasing error in the ATL08_20 product with an increasing slope for both forested and non-forested areas ( Fig. 8c and d).
Over flat terrain (with slope lower than 5 • ) tree height did not impact the accuracy of the ALT08_20 data (Fig. 8f). Only dense canopies (larger than 95%) influence the accuracy of the ALT08_20 product, when considering samples on flat terrain (Fig. 8h).

Accuracy of GEDI terrain height estimates
Environmental conditions (land cover, terrain slope, canopy cover) impact the accuracy and precision of GEDI terrain elevation estimates. However, the vertical errors caused by terrain slope observed here are of the same magnitude as would be expected due to horizontal geolocation error, i.e., vertical error = horizontal error*tan(slope). Due to the absence of complex vertical structure, the most accurate and precise GEDI terrain height estimates were found for non-forested areas (Table 1, Fig. 5a). While the mean error of GEDI terrain height for different temperate forest types is below 1 m, LE90 (as an indicator for precision) is much higher than that of non-forested land cover classes (Table 1). As expected, tropical forests represent the most challenging environment to estimate terrain elevation. In this environment GEDI mean error has a positive offset (i.e., overestimation of GEDI terrain elevation) with the highest mean error for tropical upland forests (>5 m). Canopies with dense cover (multiple vegetation layers or dense understory) cause a weak return from the ground. In this case, alternative algorithms with a lower signal end threshold (e.g., algorithm 5) may provide more accurate estimates than the default algorithm ( Fig. 9b and  e). However, in our tropical floodplain and tropical secondary forests, GEDI's default algorithm outperforms alternative algorithms in terms of statistical metrics (Fig. 9a and d; Fig. 9c and f). Furthermore, over floodplain forests we would recommend using only those GEDI samples that have been collected during the low flood season, since seasonal flood depths can be up to several meters. For LiDAR bathymetry, a wavelength of 1064 nm is typically used to derive water surface elevation, since at that wavelength the signal does not penetrate into the water column. A wavelength of 532 nm is used for measuring water depth from an aerial platform (Mandlburger et al., 2013), but use of this wavelength is not feasible from space (Abdallah et al., 2012). First results on the impact of flooding on terrain elevation estimation on the Juruá floodplain indicate an overestimation of elevation by around 5-10 m, when comparing GEDI elevation estimates acquired during high water to those acquired during low water stage. This seasonal variability offsets the higher accuracy expected for floodplain environments based on the flatness of the terrain. The choice of alternative algorithm settings is dependent on the vertical complexity of the study site. Nevertheless, the GEDI default algorithm works in most conditions (Table 1, Fig. 5a and b).
Our results showed that terrain slope impacts the accuracy of GEDI terrain elevation estimates for forested and non-forested areas ( Fig. 5c  and d). These errors can be partly explained by the geolocation uncertainty of GEDI at footprint level. Therefore, gridded GEDI products at 1 km are probably less affected by terrain slope. GEDI terrain deviations increase negatively, i.e., GEDI underestimates terrain elevation with larger slope probably due to multi-modal ground returns (Hancock et al., 2012). An increasing residual error of GEDI terrain estimates on steep slopes is partly caused by geolocation error of GEDI samples of around 10 m for version 2 (Beck et al., 2021). Similar trends in increased GEDI residual errors with an increasing slope were observed in temperate forests of central Germany (Adam et al., 2020), in Mediterranean temperate forests of Southwest Spain (Quirós et al., 2021), in an alpine forest region of northern Italy (Kutchartt et al., 2022), in different ecozones across the United States (Liu et al., 2021) and across Eucalyptus plantations in Brazil (Fayad et al., 2021). Kutchartt et al. (2022) reported that terrain slope is the most important factor influencing GEDI terrain height accuracy among other environmental (forest type, canopy cover) and sensor (beam sensitivity) parameters.
In contrast to Liu et al. (2021), we did not observe an impact of tall trees on GEDI terrain elevation estimates, when considering samples over flat terrain. Considering all terrain slopes, GEDI residual errors are higher for larger trees, but this is most likely driven by the slope, as the residual error did not increase with increasing vegetation height on flat terrains ( Fig. 5e and f). In contrast to canopy height, dense canopies affected GEDI terrain accuracies with the highest errors for canopies denser than 95% over flat areas. Dense canopies (both horizontally and vertically) hinder photons from reaching the ground, so that the ground signal can be weak and thus misclassified. However, canopy cover might have higher impact on residual errors depending on the definition of canopy cover, e.g., if canopy cover is defined as percentage of returns above 2 m rather than 5 m.

Accuracy of ICESat-2 terrain height estimates
We observed similar trends for ICESat-2 ATL08_100 and ATL08_20 terrain height estimates as for GEDI-based estimates. For instance, in forested areas ICESat-2 terrain estimates are less precise than in nonforested sites, with tropical upland forests being the most challenging environment among the forest types analyzed here. Furthermore, increasing slope leads to a negative offset of ICESat-2 terrain heights. An increase of residual errors of ICESat-2 terrain estimates was observed for both 100 m and 20 m transects. That is, the spatial size of a transect might affect the accuracy of the elevation estimate, but is probably not the main driver for this.
Canopy height itself (on flat terrain) did not impact the accuracy and precision of ICESat-2. Similarly as for GEDI-based elevation estimates, the precision of ICESat-2 terrain elevation over areas with very dense canopies (>95%) was the lowest among the canopy cover classes.
By comparing ATL08_100 and ATL08_20 products, most statistical metrics of the ATL08_20 product were more accurate across the land cover classes and forest types (Tables 1 and 2). A spatially smaller segment (e.g., 20 m) covers more homogeneous areas than a bigger segment (e.g., 100 m) leading to a reduced spatial heterogeneity and more accurate estimates. Furthermore, a higher number of samples of the ATL08_20 product reduces the impact of outliers on global statistics. Conversely, over a 100 m segment there are potentially five times more terrain photons than over a 20 m segment. However, as we can see, the ATL08_20 product provides slightly more accurate elevation estimates than the 100 m version, indicating that there are sufficient terrain photons over a 20 m transect. Unfortunately, in the current ATL08 version (v5) the number of classified terrain photons for a 20 m segment is not provided.
A first assessment of ATL08 terrain estimates version 1 in a boreal forest of Finland showed an RMSE of 0.85 m (Neuenschwander and Magruder 2019) compared to the national ALS data. Similar accuracies with RMSE of 0.73 m and ME of − 0.07 m for ATL08 terrain estimates version 3 were reported in southern Finland in . Across the six major US ecozones Malambo and Popescu (2021) estimated an ME of ATL08 terrain height product version 2 to be 0.18 m. In five of six forest types (except the tropical upland forest) we obtained ME below 1 m. Malambo and Popescu (2021) reported an optimal canopy cover for terrain height retrieval below the forest being between 40 and 70% independent of beam strength and acquisition time. In this study, however, we did not observe a difference in ICESat-2 terrain residual error for canopy densities between 0 and 80% (Figs. 7h and 8h). The reason for this is partly due to different canopy cover/density definitions, as Malambo and Popescu (2021) defined canopy cover as a percentage of returns above 2 m rather than 5 m.

Data screening and inter-comparison of sensor performance
The accuracy of ICESat-2 and GEDI terrain height estimates can be further improved by applying additional information available in the datasets (e.g., number of modes in GEDI, number of terrain photons and estimated terrain uncertainty in ICESat-2), setting maximum deviation between estimated terrain elevation and global DEMs (e.g., 30 or 50 m), although the number of samples is thereby reduced. Furthermore, to increase precision of GEDI terrain estimates further filters can be applied, e.g., removing samples where the elevation difference among the six algorithm setting groups is higher than a certain threshold. For instance, Potapov et al. (2021) applied a threshold of 2 m among the six algorithms. However, this can bias the samples to the areas of lower canopy cover, and thus, bias average properties of an area. Data screening based on beam strength (strong-power/weak-coverage) and time of acquisitions (day/night) are important for canopy height retrieval (Liu et al., 2021;Malambo and Popescu 2021), and thus, for aboveground biomass estimation (e.g., Duncanson et al., 2020), but impacts to a lower extent the estimation of terrain elevation (Liu et al., 2021;Tian and Shan 2021). Local environmental conditions should be considered when using spaceborne LiDAR-based terrain estimates, e.g., selecting data acquired during low-water season over floodplain forest areas. Furthermore, the impact of leaf-on/-off conditions can be investigated for both terrain and canopy height estimates. Finally, current versions of both ICESat-2 and GEDI terrain elevation estimates are impacted by terrain slope and to a lower extent by canopy cover. Thus, ICESat-2 and GEDI samples collected over flat terrain are more reliable.
The objective of the study is not to compare performance between the spaceborne sensors, but rather to provide statistics of current available versions under different conditions (land cover, canopy cover, terrain slope). The goodness-of-fit statistics are not directly comparable between the sensors, since they possess a different footprint size. The difference in scale might be a strong driver for statistics. For instance, variability of topography and of canopy height and density are likely to be greater over a 100 × 14 m segment than over a 25 m circle. Furthermore, an estimate for an area of 1400 m 2 will naturally differ from an estimate for an area of 280 m 2 (20 × 14 m) or ~491 m 2 (GEDI footprint). One can assess the effects of scale by e.g., averaging three GEDI footprints (sampled areã1473 m 2 ), which are comparable to the area of 1400 m 2 sampled by ICESat-2.
One important point in assessing spaceborne-based terrain elevation estimates is the quality of the reference data. In the temperate forests and in savannas the accuracy of reference DTMs is presumably high, owing to high point cloud density, rather homogeneous even-aged forest stands with few tree species, or single scattered trees in savanna. In the tropical forest, where the ALS point density was lower and the forests are more complex with denser canopies and multiple vegetation layers, the accuracy of the ground return might be lower than in temperate forests. Both spaceborne LiDAR products are still in their early versions and it is expected that they will be continuously improved in terms of terrain elevation estimates.

Conclusions
ICESat-2 and GEDI data provide freely available near-global threedimensional data of the Earth's surface. Here, we examined the accuracy of the terrain elevation estimates of the earlier versions of both LiDAR datasets across four land cover and six forest types using high-resolution airborne data as a reference. Specifically, we analyzed the accuracy of GEDI L2A version 2 and the ICESat-2 ATL08 version 5 product for 100 m and 20 m segments in different ecosystems and assessed the impacts of environmental conditions such as terrain slope, canopy height and canopy cover, on their accuracies. Our results indicated that both LiDAR missions provided accurate terrain elevation estimates across different land cover classes and forest types with ME lower than 1 m, except for tropical forests. Here, however, using GEDI alternative algorithms could improve accuracy of terrain elevation estimates. Specific environmental (e.g., terrain slope, phenological state, canopy cover) and sensor (e.g., GEDI terrain estimation algorithm, GEDI beam sensitivity, ICESat-2 number of terrain photons) parameters can be analyzed to identify more confident ICESat-2 and GEDI-based terrain estimates.
Our results showed that terrain slope negatively impacts the accuracy of both ICESat-2 and GEDI terrain elevation estimates for forested and non-forested areas. In case of GEDI the impact of slope is, however, partly attributed to the current horizontal geolocation accuracy (10 m). The vertical errors caused by terrain slope observed here are of the same magnitude as would be expected due to horizontal geolocation error (i.e. vertical error = horizontal error*tan(slope)). Further, dense canopies affected the accuracy of spaceborne LiDAR terrain estimates, while canopy height itself did not. In future versions of the dataset the horizontal geolocation accuracy will be improved, so that accuracy of terrain estimates will presumably be enhanced. Overall, both datasets provide new direct measurements of vertical vegetation structure in undersampled regions such as tropical forests and Siberian boreal forests. These datasets can be further fused with optical (Healey et al., 2020;Narine et al., 2019;Potapov et al., 2021) and radar (Qi et al., 2019) imagery to provide wall-to-wall estimates of forest height and aboveground biomass and to produce regional bare-earth DTMs.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.