New ICESat-2 Satellite LiDAR Data Allow First Global Lowland DTM Suitable for Accurate Coastal Flood Risk Assessment

: No accurate global lowland digital terrain model (DTM) exists to date that allows reliable quantiﬁcation of coastal lowland ﬂood risk, currently and with sea-level rise. We created the ﬁrst global coastal lowland DTM that is derived from satellite LiDAR data. The global LiDAR lowland DTM (GLL_DTM_v1) at 0.05-degree resolution (~5 × 5 km) is created from ICESat-2 data collected between 14 October 2018 and 13 May 2020. It is accurate within 0.5 m for 83.4% of land area below 10 m above mean sea level ( + MSL), with a root-mean-square error (RMSE) value of 0.54 m, compared to three local area DTMs for three major lowland areas: the Everglades, the Netherlands, and the Mekong Delta. This accuracy is far higher than that of four existing global digital elevation models (GDEMs), which are derived from satellite radar data, namely, SRTM90, MERIT, CoastalDEM, and TanDEM-X, that we ﬁnd to be accurate within 0.5 m for 21.1%, 12.9%, 18.3%, and 37.9% of land below 10 m + MSL, respectively, with corresponding RMSE values of 2.49 m, 1.88 m, 1.54 m, and 1.59 m. Globally, we ﬁnd 3.23, 2.12, and 1.05 million km 2 of land below 10, 5, and 2 m + MSL. The 0.93 million km 2 of land below 2 m + MSL identiﬁed between 60N and 56S is three times the area indicated by SRTM90 that is currently the GDEM most used in ﬂood risk assessments, conﬁrming that studies to date are likely to have underestimated areas at risk of ﬂooding. Moreover, the new dataset reveals extensive forested land areas below 2 m + MSL in Papua and the Amazon Delta that are largely undetected by existing GDEMs. We conclude that the recent availability of satellite LiDAR data presents a major and much-needed step forward for studies and policies requiring accurate elevation models. GLL_DTM_v1 is available in the public domain, and the resolution will be increased in later versions as more satellite LiDAR data become available. of


Introduction
Coastal flood risk assessments typically require a DTM accuracy of at least 0.5 m [1], given that the range of tidal fluctuations and sea-level rise is measured in a few meters only [2]. The uncertainties resulting from lower DTM accuracy have profound impacts on the usability of coastal land suitability assessments and on the confidence in sea-level rise impact projections, amongst others.
At present, DTMs that present the true bare ground surface with such accuracy exist for only a few coastal lowland regions, mostly in affluent countries because of the high cost involved. For most of the world's coastal lowlands, global digital elevation models (GDEMs) are used instead, with SRTM being the oldest and most commonly applied [3][4][5][6][7]. The vertical accuracy of SRTM is, however, very low at 16 m (90% confidence level; RMSE = 9.73 m; [8]) at the global level and higher but

Generating a Global Grid and Land Mask at 0.05-Degree Resolution
A global grid of 0.05 × 0.05 degrees is created that measures~5.6 × 5.6 km near the equator and~2.8 × 5.6 km at 60 degrees latitude, i.e., roughly 5 × 5 km overall. We use version 3.6 of the global administrative vector file available from the gadm database (www.gadm.org) to create a land mask. A grid cell that contained more than 50% 'land' is fully considered land. Furthermore, only when cells consisted of water for less than 50% according to the SRTM water body dataset [32], they are considered land. All maps used and presented in this study are resampled to this grid and corrected for this land mask before further analysis.

ICESat-2 Data
Version 3 of the ICESat-2 ATL08 geophysical data product over the period 14 October 2018-13 May 2020 is obtained from the National Snow and Ice Data Center [33]. The algorithm producing the ATL08 product extracts terrain and canopy heights from the ATLAS point clouds [27]. Photons are classified as ground, canopy, or noise. After classification, parameters, such as terrain and canopy height, are calculated along fixed 100 m data segments. A minimum of 50 signal photons (i.e., photons reflected off a surface) are required within the 100 m data segment for the algorithm to return the terrain parameters [27]. We have extracted the median terrain height (h_te_median; median of the photon heights above the WGS84 ellipsoid, classified as terrain within the segment) at 100 m segments from each of the 88,127 currently available ground tracks.

Local DTMs
To determine the accuracy of the GLL_DTM_v1 created in this study, existing well-described and accurate DTMs for major lowland regions across three continents are used: (1) The Everglades in the USA, (2) The Netherlands, and (3)   For the Everglades, we use the NOAA sea-level rise DEM, a LiDAR-derived DEM for coastal areas at 5 m resolution, which is available from [34].
For the Netherlands, we use the national DTM (Actueel Hoogtebestand Nederland), AHN3 [35] that is created from LiDAR data at a 5 m resolution, which is available from [36].
The Mekong Delta dataset [19] is a non-LiDAR DTM that is generated from a topographic map at 500 m resolution, using nearly 20,000 ground elevation measurements, which is available from [37].

Datum Conversion and Referencing to MSL
Both ICESat-2 and TanDEM-X data are referenced to the WGS84 ellipsoid, while SRTM, MERIT, and CoastalDEM elevations are orthometric heights referenced to the EGM96 geoid. For comparison, ICESat-2 and TanDEM-X elevations are transformed to the EGM96 geoid using the 5-min version available from https://sourceforge.net/projects/geographiclib/files/geoids-distrib/.
Datum levels for the local reference DTMs are also transformed. For the Everglades dataset, the GEOID12B [44] correction is used in the VDATUM software, version 4.0.1 [45]. For the AHN3, the RDNAPTRANS 2018 parameters, in combination with PROJ6 [46], are used. Following these local transformations to the WGS84 ellipsoid, the DTMs are referenced to the EGM96 geoid.
To convert the vertical datum from the EGM96 geoid to MSL, we use the mean dynamic topography (MDT), which is the difference between the mean sea surface and the geoid. We use an estimate by [47], who calculated the MDT by combining geodetic data (i.e., altimetric mean sea surface over the period 1993-2012 and an accurate geoid) with in situ data, at 0.25-degree resolution (http://www.aviso.altimetry.fr/). For referencing to MDT on land, the closest value at sea is used. The Mekong Delta dataset is already referenced to local MSL, so no further conversion is necessary.

Generating GLL_DTM_v1
We observe that some of the ICESat-2 ground tracks contain erroneous extreme outlier values. Since the ATL08 algorithm is still in development, and such values are likely to be removed in the future [27], we do not attempt to correct for these erroneous data, but rather apply a hard cut-off value to exclude them. All data segments with a median terrain height below −7 m +MSL are removed because they are considered unrealistic in coastal lowland, noting that −7 m +MSL is the lowest land elevation value found in pumped polders in the Netherlands. Visual inspection confirms that nearly all erroneous data segments are indeed removed in this way.
Since ICESat-2 data coverage is insufficient to create a DTM of similar resolution as the GDEMs, we resample the 100 m data segment elevation values to a 0.05-degree land grid by taking the median of 100 m data segments within each grid cell. The 0.05-degree grid resolution is found to be the most suitable resolution at the moment, given that nearly all (98.8%) of such cells have terrain elevation data from at least one ICESat-2 data segment. A total of 97.6% of 0.05-degree grid cells have more than five data segments. The distribution of the number of data segments within each 0.05-degree lowland cell is provided in Figure 3. To create an approximate full coverage lowland DTM, values in 'no data' cells (1.2% of the cell total) are interpolated between data cell values through inverse distance interpolation. This interpolation method is applicable in areas where landscape features are gradual, as is usually the case in lowlands that tend to be very flat [25]. For each region in the world, excluding Antarctica ( Figure 2), areas within specified elevation classes are calculated applying the equal area projection [48,49] and compared with the different GDEMs.

Accuracy Assessment by Comparison with Local DTMs
After transformations, both the local DTMs and GDEMs are resampled to the same 0.05-degree grid as the ICESat-2 data, by determining the median value within cells. The difference between GLL_DTM_v1 and the local DTMs at 0.05-degree resolution is mapped for each area, and the mean error (Equation (1)), mean absolute error (Equation (2)), and RMSE (Equation (3)) are calculated and tabulated. This is done separately for land, that is, below 2, 5, and 10 m +MSL according to the local DTMs. Similarly, also the GDEMs are compared with the local DTMs at 0.05-degree resolution.
where N is the total number of samples, ∆h i = elevation difference between assessed DTM and reference DTM at sample point (i).

GLL_DTM_v1 Accuracy Compared to GDEMs
The accuracy assessment covers 74,604 km 2 of the lowland area (<10 m +MSL). The GLL_DTM_v1 is found to be accurate within 1 m for 95.5% and within 0.5 m for 83.4% of land area below 10 m +MSL, with a mean error of 0.06 m, a mean absolute error of 0.32 m, and an RMSE of 0.54 m, when compared with three local DTMs (Table 1, Figure 4). The RMSE is similar for the three comparison areas (range 0.44-0.67 m), indicating a high degree of consistency across regions in the GLL_DTM_v1 product. Table 1. Accuracy of GLL_DTM_v1 and existing GDEMs for areas lower than 2, 5, and 10 m +MSL (as compared to local DTMs). Presented are mean elevation differences between the local DTM, GLL_DTM_v1, and existing GDEMs, as well as RMSE and the area percentage differences within ranges −1-+1 m and −0.5-+0.5 m.    Figure 5, Figures S1 and S2). Mean errors are similar to mean absolute errors across GLL_DTM_v1 and GDEMs, except for CoastalDEM, where the mean error averaged over the three areas is negative at −0.32 m.  Figures S1 and S2, respectively.
For areas below 2 m +MSL, the difference between GLL_DTM_v1 and existing GDEMs is greatest, with extent according to GLL_DTM_v1 being 295% of that of SRTM90, 208% of MERIT, 111% of CoastalDEM, and 158% of TanDEM-X (Table 2, Figure 6). For global land below 10 and 5 m +MSL, the GLL_DTM_v1 extent is also greater at 151-201%, 105-121%, 108-109%, and 117-132%, respectively (Table 2, Figure 6). Table 2. Coastal lowland areas (in 10 3 km 2 ) below 2, 5, and 10 m +MSL in different regions (60N-56S, as covered by SRTM), applying the GLL_DTM_v1 generated from ICESat-2 data in this study and other GDEMs. The region definition is shown in Figure 2.   The coastal lowland extent, as determined from GLL_DTM_v1, is largest in the greater Southeast Asia region, at 36%, 36%, and 35% of the global total at elevations below 10, 5, and 2 m +MSL. For land below 2 m +MSL, this is followed by South America at 22%, Africa and the Middle East at 15%, North America at 12%, Europe and Russia at 6%, East and Central Asia at 5%, and Australia and Oceania at 4% ( Table 2).
The variation between sources in lowland area below 2 m +MSL is greatest in South America, with the extent of SRTM90, MERIT, CoastalDEM, and TanDEM-X being 14%, 29%, 71%, and 39% that of GLL_DTM_v1. For Europe and Russia, these numbers are 101%, 85%, 129%, and 125%, presenting the smallest variation between sources ( Table 2).
The lowland area in all of Asia, i.e., along the coastlines of greater Southeast Asia and East and Central Asia combined, is 44% of global extent below 10 m and 5 m +MSL and 40% of extent below 2 m +MSL, according to GLL_DTM_v1.

Impact of DTM Accuracy on Lowland Extent Estimation
The low errors and high consistency across regions in GLL_DTM_v1 compared to radar-based GDEMs demonstrates that the ICESat-2 satellite LiDAR data provide a unique opportunity to improve global lowland area mapping with substantially higher accuracy than has been possible so far and, therefore, will allow more accurate environmental and climate impact assessments relying on such data products.
The higher global extent of coastal lowland below 10, 5, and 2 m +MSL, determined by GLL_DTM_v1 compared to other GDEMs ( Table 2), confirms that earlier radar-based GDEMs have underestimated lowland area because of higher apparent land surface elevation, caused at least in part by radar being unable to fully penetrate the vegetation canopy [8,16,17]. The difference between the sources is greatest at the lowest elevation range, below 2 m +MSL, because there the relative effect of vegetation and referencing inaccuracies in GDEMs may be most pronounced.
Of the four radar-based GDEMs assessed, TanDEM-X has the highest accuracy, as has also been found in other recent studies comparing TanDEM-X with SRTM and MERIT [17,50], but the error within 0.5 m over 37.9% of extent is still far below that of GLL_DTM_v1 at 83.4%.
The advantage of using satellite LiDAR data for creating DTMs is especially clear when considering densely vegetated areas. The lowlands of Papua and West Papua Provinces in Indonesia, for instance, that are still largely covered with forest, cover 29,180 km 2 of land below 2 +MSL according to GLL_DTM_v1 (Figure 7), but only 10.0% of that (2911 km 2 ) according to TanDEM-X, and 0.7-72.8% (215-21,233 km 2 ) according to SRTM90 and derived products MERIT and CoastalDEM. Similarly, for the Amazon Delta, GLL_DTM_v1 finds 61,217 km 2 of land below 2 m +MSL, but only 19.0% of that (11,626 km 2 ) according to TanDEM-X, and 2.8-44.3% (1723-27,097 km 2 ) according to the other three GDEMs. Where the radar-based GDEMs detect lowland in these areas, it is generally in gaps in the forest cover, as is clear in Figure 7. Of the GDEMs, MERIT best approaches the shape of lowland areas presented by GLL_DTM_v1 in Papua and Amazon Delta, but it presents a land surface that is several meters too high.

Effect of Resolution on the Accuracy
To test the effect on the accuracy of aggregating data at 0.05-degree spatial resolution, we compare the result with that at the native (original to source) resolution. The distribution of differences between the GDEMs compared in this study and the three local DTMs at both native and 0.05-degree resolutions are shown in Figure 8. Similar to GLL_DTM_v1, the original ICESat-2 segment data have far higher accuracy than the other GDEMs compared to the three local DTMs with average differences ranging between −0.11 and 0.10 m and RMSE between 0.92 and 1.17 m (Table 3). These RMSE values are similar to the RMSE of 0.85 m found by Neuenschwander and Magruder [53] for ATL08 terrain heights when comparing an ICESat-2 track over vegetation in Finland with spatially coincident airborne LiDAR data.  Table 3. Accuracy of ICESat-2 data segments and existing GDEMs at native (90-100 m) spatial resolution for areas lower than 10 m +MSL (as determined from the local DTMs). Presented are mean elevation differences between the local DTM and ICESat-2 data segments and existing GDEMs as well as RMSE. Differences for SRTM90, MERIT, CoastalDEM, and TanDEM-X are computed only when coinciding with ICESat-2 data segments. The tails of the boxen plots ( Figure 8) show that outlier values occur in the original ICESat-2 data segments. This illustrates that direct interpolation between individual segments to generate a DTM would result in erroneous DTM cell values and that the impact of outlier values on the DTM is reduced by aggregating to larger cells.
From Figure 8, it is also clear that the GDEMs at native resolution all display much longer tails, i.e., contain more outlier values, than the ICESat-2 data. This confirms that the ICESat-2 LiDAR data is intrinsically more accurate than radar-based sources. Aggregated to 0.05-degree resolution in GLL_DTM_v1, the range of ICESat-2-based DTM difference values is several times smaller than that of any GDEM for any of the three local comparison areas.

Impact of DTM Accuracy on Flood Risk Assessments
An indication of the impact of differences between DTMs and GDEMs on assessing flood risk is obtained by demonstrating the depth of inundation that would occur if water levels are 2 m above current MSL (Figure 9). Where the local TOPODEM DTM and GLL_DTM_v1 yield limited land area with a flood depth more than 2 m, at 0.5% and 4.7% of the Mekong Delta extent, respectively, this is 24.8%, 1.8%, 79.1%, and 0.3% for SRTM90, MERIT, CoastalDEM, and TanDEM-X. Such large differences are meaningful to flood risk assessments as risk is not only determined from flood extent but also the depth of inundation. Flood depths of 2 m or more are generally considered life-threatening, whereas depths of 1 m or less will mostly damage infrastructure and crops.

Conclusions
We show that satellite LiDAR data can yield a global lowland DTM that closely matches three local DTMs, across three continents, within 0.5 m over 83.4% of their area lower than 10 m +MSL, and far exceeds the accuracy of existing radar-based GDEMs below 10 m +MSL that are accurate within 0.5 m for 37.9% (TanDEM-X) at best. This confirms that satellite LiDAR data presents a way forward towards creating accurate global lowland DTMs that will be suitable for assessments of current and future flood risk and land management options.
We demonstrate that the availability of a more accurate global lowland DTM will have major implications for flood risk assessments, in terms of the extent and nature of land at the highest risk. The global land area below 2 m +MSL, often defined as land that is already inundated periodically, even without sea-level rise, is around 1 million square kilometers, which is equivalent to 25 times the land area of the Netherlands. This is more than is calculated using any other source and almost three times more than what is found using SRTM90 data that are used in most flood risk assessments to date.
The greatest underestimation of lowland extent by existing radar-based GDEMs occurs in remaining forested areas, such as Papua and Amazon Delta. With GLL_DTM_v1, it is possible to more accurately determine the societal, economic, and ecological risks of forest clearing for the development of agriculture or urban expansion in such areas.
The result presented here applies to a relatively coarse 0.05-degree (~5 × 5 km) resolution DTM, based on 19 months of ICESat-2 data collection. As additional satellite LiDAR will become available, from ICESat-2 and other sources, such as GEDI [55], it will be possible to further improve both the resolution and accuracy of DTMs, increasing usability in detailed assessments. However, for the application of identifying large lowland areas at the regional and global scale, the current resolution is appropriate.
We recommend a shift from radar-based DEMs to LiDAR-based DTMs for lowland areas. The GLL_DTM_v1 product can be downloaded from https://doi.org/10.17632/v5x4vpnzds.1. GLL_DTM updates at a higher resolution will be uploaded at this location.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/17/2827/s1, Figure S1: Comparison of GLL_DTM_v1 and GDEMs with the local DTM for the Everglades, Figure S2: Comparison of GLL_DTM_V1 and GDEMs with the local DTM for the Netherlands.  Acknowledgments: The idea to create GLL_DTM_v1 originated in discussions with Marcus Wishart, George Stirrett Wood, and Abedalrazq Khalil at the World Bank. We acknowledge the use of ICESat-2 data from the National Snow and Ice Data Center. Bregje van Wesenbeeck is gratefully acknowledged for constructive comments on an earlier version of the manuscript. We wish to thank the anonymous reviewers for their constructive and useful comments and suggestions that enabled us to improve the content of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: