Perspectives on Digital Elevation Model (DEM) Simulation for Flood Modeling in the Absence of a High-Accuracy Open Access Global DEM

Open-access global Digital Elevation Models (DEM) have been crucial in enabling ﬂood studies in data-sparse areas. Poor resolution ( > 30 m), signiﬁcant vertical errors and the fact that these DEMs are over a decade old continue to hamper our ability to accurately estimate ﬂood hazard. The limited availability of high-accuracy DEMs dictate that dated open-access global DEMs are still used extensively in ﬂood models, particularly in data-sparse areas. Nevertheless, high-accuracy DEMs have been found to give better ﬂood estimations, and thus can be considered a ‘must-have’ for any ﬂood model. A high-accuracy open-access global DEM is not imminent, meaning that editing or stochastic simulation of existing DEM data will remain the primary means of improving ﬂood simulation. This article provides an overview of errors in some of the most widely used DEM data sets, along with the current advances in reducing them via the creation of new DEMs, editing DEMs and stochastic simulation of DEMs. We focus on a geostatistical approach to stochastically simulate ﬂoodplain DEMs from several open-access global DEMs based on the spatial error structure. This DEM simulation approach enables an ensemble of plausible DEMs to be created, thus avoiding the spurious precision of using a single DEM and enabling the generation of probabilistic ﬂood maps. Despite this encouraging step, an imprecise and outdated global DEM is still being used to simulate elevation. To fundamentally improve ﬂood estimations, particularly in rapidly changing developing regions, a high-accuracy open-access global DEM is urgently needed, which in turn can be used in DEM simulation.


INTRODUCTION
Digital Elevation Models (DEM) are a gridded digital representation of terrain, with each pixel value corresponding to a height above a datum.Since the pioneering work of Miller and Laflamme (1958), DEMs have grown to become an integral part of a number of scientific applications.DEMs can be created from ground surveys, digitizing existing hardcopy topographic maps or by remote sensing techniques.DEM's are now predominantly created using remote sensing techniques with Smith and Clark (2005) observing the benefits that a large spatial area can be mapped by fewer people at a lower cost.Remotely sensing techniques include photogrammetry (Uysal et al., 2015;Coveney and Roberts, 2017), airborne and spaceborne Interferometric Synthetic Aperture Radar (InSAR) and Light Detection And Ranging (LiDAR).Spaceborne InSAR is the most common technique to create global DEMs and is the technology behind the most widely used open-access global DEM; the Shuttle Radar Topography Mission (SRTM).An overview of free and commercial global DEMs is given in Figure 1.Despite being acquired in 2000, SRTM is still the most popular global DEM because of its accessibility, feature resolution, vertical accuracy and a lower amount of artifacts and noise compared to alternative global DEMs (Rexer and Hirt, 2014;Jarihani et al., 2015;Sampson et al., 2016;Hu et al., 2017).Yet, recently released products such as MERIT and TanDEM-X 90 could change that.Commercial DEMs (Planet Observer, 2017;Takaku and Tadono, 2017;InterMap, 2018) are often prohibitively costly and have a lack of comparison studies.Airborne LiDAR has a higher precision and accuracy owing to its ability to penetrate vegetation and its reduced vulnerability to scatter but is largely limited to a handful of countries and can be expensive to acquire.These characteristics are conducive in creating a high-quality (vertical error <1 m) 'bare-earth DEM, ' where objects (e.g., buildings and vegetation) have been removed from the elevation model.Such bare-earth DEMs are essential for applications, such as flood modeling, that rely on the accurate derivation of surface characteristics (e.g., slope).

Characteristics of DEM Errors
DEM errors occur in both the horizontal and vertical directions.Errors propagate from the input data used in creating a DEM right through to calculating surface derivatives and using DEMs in complex applications (Hutchinson and Gallant, 2000;Fisher and Tate, 2006).Whatever their source, DEMs can appear to provide a definitive and plausible representation of topography which can often lull the user into a false sense of security regarding their accuracy, with many users unaware of DEM errors or how to treat them (Wechsler, 2003(Wechsler, , 2007)).Whilst accuracy statistics such as RMSE provide an indication of DEM accuracy, they assume error to aspatial (Hunter and Goodchild, 1997;Carlisle, 2005;Fisher and Tate, 2006;Wechsler, 2007).Invoking Tobler's First Law of Geography, whereby he noted that "nearby things are more similar than distant things" (Tobler, 1970), we know error varies spatially so DEM error is spatially autocorrelated.Indeed, Holmes et al. (2000) observe that "although global (average) error is small, local error values can be large, and also spatially correlated."The spatial variation of DEM error is most frequently estimated by calculating accuracy statistics of areas disaggregated by slope and/or landcover class, and more rarely spatial structure of error.Wise (2000) categorized DEM errors as systematic, blunders or random.These types of errors derive from: (a) deficient spatial sampling and/or age of data; (b) processing errors such as interpolation or numerical errors; (c) measurement errors from poor positional inaccuracy, faulty equipment or observer bias (Wechsler, 2007).Systematic errors occur in the DEM creation procedure by processing techniques that can cause bias or artifacts.Blunders arise from human error (Wise, 2000) or equipment failure (Fisher and Tate, 2006).Random errors occur in any system of measurement due to the wealth of measurement and operational tasks performed to create a DEM (Wise, 2000;Fisher and Tate, 2006), and remain even after known blunders and systematic errors are removed (Wechsler, 2007).An example of random error is speckle noise (multiplicative noise in a granular pattern) (Rodriguez et al., 2006;Farr et al., 2007).Sources of systematic errors and blunders relevant to flood modeling derive from interpolation techniques (Desmet, 1997;Wise, 2007;Bater and Coops, 2009;Guo et al., 2010), erroneous sink filling (Burrough and McDonnell, 1998), hydrological correction (Callow et al., 2007;Woodrow et al., 2016), deficient spatial sampling causing urban features not to be resolved (Gamba et al., 2002;Farr et al., 2007), slope and aspect (foreslope vs backslope) (Toutin, 2002;Falorni et al., 2005;Shortridge and Messina, 2011;Szabó et al., 2015), striping caused by instrument setup (Walker et al., 2007;Tarakegn and Sayama, 2013) and vegetation (Carabajal and Harding, 2006;Hofton et al., 2006;Shortridge, 2006;Weydahl et al., 2007;LaLonde et al., 2010).
The aforementioned errors propagate into errors in surface derivatives including, but not limited to, slope (Holmes et al., 2000), aspect (Januchowski et al., 2010), curvature (Wise, 2011), drainage basin delineation (Oksanen and Sarjakoski, 2005) and upslope contributing area (Wu et al., 2008).As many models rely on these surface derivatives (e.g., change in slope is the dominant control on flow in flood models), error propagation from DEMs can substantially affect results of models that use these surface derivatives.Yet, taking flood models as an example, sensitivity analysis has largely focussed on hydraulic parameters and has under-represented DEM errors (Wechsler, 2007).Davis and Keller (1997) aptly sum up the problem of DEM error with their remark that 'landscapes are not uncertain, but knowledge about them is.'

Flood Inundation Models and DEM Error
Topography is arguably the key factor for the estimation of flood extent (Horritt and Bates, 2002), but typically flood models use a limited number of DEMs and instead choose to explore the uncertainty associated with other hydraulic parameters (Wechsler, 2007).Studies that do use multiple DEMs either resample DEMs to a coarser resolution to explore the effect of resampling strategies and/or scale (Horritt and Bates, 2001;Neal et al., 2009;Fewtrell et al., 2011;Saksena and Merwade, 2015;Savage et al., 2016b;Komi et al., 2017), or compare flood extents using different DEM products (Li and Wong, 2010;Jarihani et al., 2015;Bhuyian and Kalyanapu, 2018).Generally speaking, the quality of flood predictions increases with higher resolution DEMs.Higher resolution DEMs are more important when modeling urban environments (Fewtrell et al., 2008) so buildings can be captured.Resolution can be less important for rural environments with Savage et al. (2016a) concluding that running simulations finer than 50 m had little performance gain without occurring additional unnecessary computational cost.Too much detail can induce spuriously precise results which does not represent the uncertainties in making flood predictions (Dottori et al., 2013;Savage et al., 2016a).In data-sparse regions, a limited number of global DEM products dictates that only a single DEM is used, with this most commonly being SRTM (Yan et al., 2015).Whilst understandable, using a single DEM leads to a dangerous situation where spuriously precise estimates of flood extent are presented which do not assess the impact of uncertain topography.DEM simulation overcomes this obstacle by making available a catalog of statistically plausible DEMs at the native resolution of the global DEM from which it is simulated from.

CURRENT ADVANCES -CORRECTING DEM ERROR
Here we identify three categories of approaches to correct DEM error: (1) DEM editing; (2) New DEMs created with improved sensing technologies and (3) Stochastic simulation of DEMs.This article focuses on the third approach.

DEM Editing
DEM error can be reduced by editing -either manually or systematically.Manual editing involves changing pixel values based on additional information or expert judgement.We speculate that this happens frequently but is seldom documented.Systematic editing involves applying algorithms, additional datasets and filters to reduce error.For example, a DEM can be hydrologically corrected through algorithms such as AGREE (Hellweger, 1997), ANUDEM (Hutchinson, 1989), outlet breaching (Martz and Garbrecht, 1999), Priority-Flood (Barnes et al., 2014) and stream burning (Saunders, 1999).Namely, the HydroSHEDs global hydrography dataset makes use of hydrological correction techniques to create invaluable maps such as flow direction, river networks and catchment masks (Lehner et al., 2008).DEMs have also been edited to correct errors from vegetation (Baugh et al., 2013;Pinel et al., 2015;Su et al., 2015;O'Loughlin et al., 2016;Ettritch et al., 2018;Zhao et al., 2018) and to compensate for the positive bias in coastal areas due to vegetation and buildings that lead to an underestimation of coastal flood exposure [e.g., CoastalDEM (Kulp and Strauss, 2018)].
The recent release of the MERIT (Multi-Error-Removed-Improved-Terrain) DEM is the most comprehensive error removal from SRTM to date.Errors are reduced by separating and removing absolute bias, stripe noise, speckle noise and vegetation bias, with the most significant improvements reported in flat regions (Yamazaki et al., 2017).Compared to SRTM, MERIT has fewer artifacts (Hirt, 2018) and a better performance in flood models compared to SRTM (Chen et al., 2018).Whilst a significant improvement on SRTM, MERIT is still fundamentally based on SRTM data and is thus limited by the errors in SRTM.The next new edited DEM will be NASADEM (Crippen et al., 2016) set to be released in 2018, but again this is a reprocessing of SRTM.

New DEMs
The most widely used open-access global DEM (SRTM) is almost two decades old.Advances in satellite technology, image processing and data storage capabilities, make creating new, more accurate DEMs entirely possible.Two new global DEMs have been recently released, namely optically derived ALOS AW3D30 (Tadono et al., 2014) (∼30 m resolution) and the SAR derived TanDEM-X 90 (Rizzoli et al., 2017) (∼90 m resolution).Both of these products are technically digital surface models, so should only be used with caution in flood models.Looking to the future, new techniques are being explored to create new DEMs.One such example of applying a new technique is the creation of a 2 m pan-Arctic DEM (ArcticDEM)1 using stereo auto-correlation techniques to overlap pairs of high-resolution optical imagery.Additionally, Ghuffar (2018) demonstrated that a 5 m DEM can be generated from Planet Labs cubesat derived PlanetScope imagery using Semi Global Matching.Alternatively, existing DEMs can be fused together to create new products (de Ferranti, 2014;Yue et al., 2017;Pham et al., 2018).For example ASTER and SRTM have been fused together to create the global EarthEnv DEM (Robinson et al., 2014).High Resolution (<10 m) open-access LiDAR data is becoming increasingly available through initiatives such as OpenTopography2 , with New Zealand the latest country to release LiDAR data for free.Despite this encouraging step, we optimistically estimate openaccess LiDAR data covers just 0.005% of the earth s land area based on data from OpenTopography and an extensive search of national mapping agencies.Global LiDAR coverage is some way off, with the limited amount of LiDAR data that is currently available almost exclusively found in developed countries.

DEM Simulation
Stochastic simulation assumes that a DEM is only a single realization amongst a host of potential realizations.DEMs are simulated by altering pixel values in accordance with the spatial error structure.A single true DEM is not created.Instead realizations provide a bound within which the true value is likely to lie.Therefore, DEM error is not reduced as such, but the bounds of error are identified.This idea is relatively well known in the field of geostatistics (Goovaerts, 1997;Hunter and Goodchild, 1997;Deutsch and Journel, 1998;Kydriakidis et al., 1999;Holmes et al., 2000).Using an ensemble of simulated DEMs has been shown to greatly affect the characterization of surface derivatives (Fisher, 1991;Veregin, 1997;Holmes et al., 2000;Endreny and Wood, 2001;Raaflaub and Collins, 2006), landslide hazard (Davis and Keller, 1997;Murillo and Hunter, 1997;Darnell et al., 2008) and flood inundation estimation (Wilson and Atkinson, 2005;Hawker et al., 2018).Research in this area, especially in the flood community, has been largely stagnant for the past decade which seems a shame given the improvement in computational resources and the number of DEMs now available.
Here we demonstrate the benefits of DEM simulation in flood inundation modeling based on the recently published work of Hawker et al. (2018).

DEM Simulation
In Hawker et al. (2018), DEM simulation is carried out by first quantifying the spatial error structure of a global DEM, and then using the fitted error covariance function to simulate plausible versions of the DEM.The fitted error covariance function was calculated for SRTM and MERIT DEMs by fitting a semivariogram to difference maps (i.e., SRTM/MERIT -reference LiDAR DEM) of 20 floodplain locations.Semi-variograms for all locations and by landcover class were produced, thus making DEM simulation possible for any floodplain SRTM/MERIT location, with the work available through the R Package DEMsimulation.For more details we refer the reader to Hawker et al. (2018).It must be noted that only a limited number of semivariograms are produced, but is nevertheless useful in creating an ensemble of DEMs.
We test the quality of DEM simulations produced by simulating MERIT DEM by landcover semi-variograms by plotting rank histograms for an ensemble of 2500 DEMs of the Ba catchment in Fiji.Rank histograms (or 'Talagrand' diagrams) (Anderson, 1996;Hamill and Colucci, 1997;Talagrand et al., 1997) are a common tool used to evaluate ensemble forecasts in meteorology, and work by ranking the verification (in our case LiDAR data) relative to the corresponding value in the ensemble in ascending order.An ideal ranked histogram is flat since the observation is indistinguishable from any ensemble member.Typically, a U-shaped rank histogram suggests undervariability in the ensemble, a dome shape over-variability, and excessive population of the extreme ranks bias.Yet, ranked histograms are notoriously difficult to evaluate and can lead to misinterpretations if done uncritically (Hamill, 2001).Nevertheless, we produce rank histograms by taking the mean of LiDAR values which fall within each ensemble pixel for the Ba catchment in Fiji (Figure 1).The rank histogram (Figure 1B) of all pixels suggests a positive bias in ensemble members as the ranks are clustered to the left.Despite the vegetation correction in MERIT, the rank histogram of mangrove covered pixels shows a large positive bias, whilst cropland has a more uniform shape.
To compensate for errors in observations (LiDAR), we added random observational noise as suggested by Hamill (2001), but this made little difference to the shape and thus is not presented here.Additionally, we compute 3 goodness-of-fit measurements: Pearson X 2 ; Jolliffe-Primo (JP) slope and JP convexity, with the null hypothesis that the rank histogram is flat (Jolliffe and Primo, 2008).These statistics confirm the stronger bias in mangroves (JP Slope) and suggest possible under-sampling with the relatively high JP convexity values.All these results are statistically significant with p-values of virtually 0.Moreover, less than 3% of pixels within the single MERIT DEM were within the error of the LiDAR (≈50 mm), whilst this was 97% for the ensembles.Therefore, the reliability of the DEM simulation is deemed satisfactory but can still suffer from systematic errors from the global DEM product being used.A higher-accuracy global DEM would therefore make this approach even more effective.

Simulated DEMs and Flood Inundation Predictions
To demonstrate the usefulness of using simulated DEMs in flood predictions we expand upon work published in Hawker et al. (2018).We simulate a total of 7500 DEMs to use in a LISFLOOD-FP Neal et al. (2012) flood model of Ba, Fiji (Figure 2A) for a 50 year return period flood event (Archer et al., 2018).The Ba catchment is predominantly agricultural floodplain, with mangroves present at the coast.DEMs are either simulated using MERIT or SRTM DEMs, and using either an average floodplain semi-variogram or semi-variograms disaggregated by landcover.
Flood predictions are compared to four models that use a single DEM -LiDAR at 30 m and 90 m resolution and MERIT and SRTM at 90 m resolution.We assume the LiDAR 30 m model is the benchmark prediction in lieu of a lack of observation data.Flood depth errors are compared against the LiDAR 30 m model are plotted for the deterministic approach using the MERIT DEM and the stochastic approach using an ensemble of simulated DEMs (Figure 2B).Whilst the DEM ensemble approach can overpredict flood extent, flood depths are often more accurate as indicated by the more neutral colors of the DEM ensemble flood map given in Figure 2B).Further analysis of predicted flood depth (Figures 2C-F) indicate the benefit of using ensembles of simulated DEMs in predicting correct water depths.For example, in location 2, the MERIT DEM does not flood, whilst the flood depth in SRTM is large (>4.8 m), but for the ensembles of DEMs the distribution of predicted flood depths are more closely aligned with the flood depths predicted in the LIDAR models.The results also highlight the differences in predictions between DEMs, so we would encourage modelers to use multiple DEMs even if DEM simulation is not used.Yet by using an ensemble of simulated DEMs, we can learn about the distribution of potential flood extent and flood depth, and thus can avoid the spurious precision when using a single DEM.

FUTURE DIRECTIONS
In this article, we have attempted to reinvigorate the idea of DEM simulation and highlight its value for flood studies.Despite repeated calls to produce a new high-accuracy open access global DEM (Schumann et al., 2014;Simpson et al., 2015), this unfortunately does not seem forthcoming.Everincreasing computing power has made even global flood simulations possible (Sampson et al., 2015), while flood modelers also often run multiple models to explore model parameter sensitivities.However, the impact of DEM error has been largely overlooked in lieu of a lack of suitable stochastic DEM data.DEM simulation overcomes this restriction, making it possible for flood modelers to use a catalog of DEMs.Working in tandem with systematic DEM editing (e.g., MERIT), DEM simulation can fill the gap until a much-needed new highaccuracy open access DEM is produced.Even when this longawaited DEM is eventually produced, DEM simulation will still be an invaluable approach for exploring the effect of DEM error in flood inundation estimates as long as good estimates of the spatial error structure can be made across a sufficient number of locations.We therefore encourage scientists to embrace geostatistics to simulate DEM ensembles and call for increased reporting of spatial dependence by DEM vendors and scientists alike.

DATA AVAILABILITY STATEMENT
The MERIT dataset can be downloaded after sending a permission request to the developer (Dai Yamazaki, yamadai@rainbow.iis.u-tokyo.ac.jp) from http://hydro.iis.utokyo.ac.jp/~yamadai/MERIT_DEM/, and is free for research and education purposes.SRTM data can be freely downloaded from https://earthexplorer.usgs.gov/.The LiDAR dataset is available from The Secretariat of the Pacific Community s Applied Geoscience and Technology Division (SPC SOPAC).LISFLOOD-FP is available for non-commercial purposes from http://www.bristol.ac.uk/geography/research/hydrology/models/ lisflood/downloads/.The code to simulate DEMs can be found at https://github.com/laurencehawker/demgenerator.

AUTHOR CONTRIBUTIONS
LH wrote the article and carried out the data analysis.PB, JN, and JR reviewed and provided comments on various drafts as well as advised on the analysis work.

FIGURE 1 |
FIGURE 1 | (A) Overview of Existing Global DEMs (free and commercial).(B) Rank Histograms for an ensemble of 2500 DEMs of the Ba catchment in Fiji, simulated from the MERIT DEM using semi-variograms of spatial error structure by landcover class.All Pixels and the two landcover classes with the most pixels (Mosaic Cropland/Natural Vegetation and Mangroves) are shown.

FIGURE 2 |
FIGURE 2 | Maximum flood water depth at four locations in Ba, Fiji for a 50-year return period event.(A) Overview of study area, with locations of four random locations to investigate differences in water depth.(B) Average flood depth for models run with MERIT DEM and simulated versions of the MERIT DEM simulated using semi-variograms per landcover class against the LIDAR 30 m model.(C) Maximum water depth distribution of each DEM ensemble simulated by different sets of semi-variograms for location 1. (D) Maximum water depth distribution of each DEM ensemble simulated by different sets of semi-variograms for location 2. (E) Maximum water depth distribution of each DEM ensemble simulated by different sets of semi-variograms for location 3. (F) Maximum water depth distribution of each DEM ensemble simulated by different sets of semi-variograms for location 4. MERIT Avg refers to MERIT DEM simulated using an 'average' floodplain semi-variogram, MERIT LC refers to MERIT DEM simulated using semi-variograms by landcover class, SRTM LC refers to SRTM DEM simulated using semi-variograms by landcover class.
LH was funded as part of the Water Informatics Science and Engineering Center for Doctoral Training (WISE CDT) under a grant from the Engineering and Physical Sciences Research Council (EPSRC), grant number EP/L016214/1.PB was supported by a Leverhulme Research Fellowship and Royal Society Wolfson Research Merit award.