Topographic uncertainty quantification for flow-like landslide models via stochastic simulations

Topography representing digital elevation models (DEMs) are essential inputs for computational models capable of simulating the run-out of flow-like landslides. Yet, DEMs are often subject to error, a fact that is mostly overlooked in landslide modeling. We address this research gap and investigate the impact of topographic uncertainty on landslide-run-out models. In particular, we will describe two different approaches to account for DEM uncertainty, namely unconditional and conditional stochastic simulation methods. We investigate and discuss their feasibility, as well as whether DEM uncertainty represented by stochastic simulations critically affects landslide run-out simulations. Based upon a historic flow-like landslide event in Hong Kong, we present a series of computational scenarios to compare both methods using our modular Python-based workflow. Our results show that DEM uncertainty can significantly affect simulation-based landslide run-out analyses, depending on how well the underlying flow path is captured by the DEM, as well as further topographic characteristics and the DEM error's variability. We further find that in the absence of systematic bias in the DEM, a performant root mean square error based unconditional stochastic simulation yields similar results than a computationally intensive conditional stochastic simulation that takes actual DEM error values at reference locations into account. In all other cases the unconditional stochastic simulation overestimates the variability of the DEM error, which leads to an increase of the potential hazard area as well as extreme values of dynamic flow properties.


Introduction
Landslides are natural hazards that occur frequently all around the world causing casualties, economic devastation, and environmental destruction. Most often, they are naturally driven, e.g. by means of long-lasting and/or intensive precipitation events, or induced by earthquakes. Yet, landslides might also be triggered or its susceptibility increased as a result of human activities, e.g. deforestation and construction. According to the United Nations Office for Disaster Risk Reduction and the Center for Research on the Epidemiology of Disasters, 378 recorded landslides from 1998 to 2017 affected 4.8 million people and caused 18414 deaths as well as several billion US dollars of economic losses [1]. Froude and Petley [2] reported that in total 55997 people were killed during 4862 fatal non-seismic landslide events from January 2004 to December 2016. Still, it has to be assumed that the damage potential of landslides is underestimated as 1) events have been underreported for decades, especially in developing countries, and 2) losses caused by co-seismic landslide events tend to be classified as secondary losses due to earthquakes. Rapid flow-like landslides, such as rock avalanches and debris flows, show a particularly high hazard potential due to their high mobility, long travel distance and fast propagation speed. In recent years, the geo-hazard community put a lot of effort into developing computational run-out models in order to assess and predict risks associated with rapid landslides and to develop mitigation strategies. Most of the models in practical use are based on a (computationally efficient) 'shallow flow type' process description and depth-averaging techniques, e.g. [3,4,5,6,7]. In these, the flowing material is treated as an 'equivalent fluid' and governed by idealized internal and basal rheologies [3]. Alternative (computationally demanding) models aim at a direct description of fully three-dimensional flow behavior. They hence offer a higher process complexity level, e.g. [8,9], yet are typically not feasible for practical hazard mitigation purposes. Detailed reviews of computational run-out models for rapid, flow-like landslide models have been published by McDougall [10] and Pastor et al. [11]. An indispensable input to any of these computational landslide run-out models is data that represents the terrain in which the slide is likely to occur. Pioneered by Miller and Laflamme [12], digital elevation models (DEMs) have become the most popular form of representing topographies in the scientific community. Methods for generating DEMs have evolved rapidly over decades from conventional approaches like field surveying and topographic map digitizing, to passive and active remote sensing, such as stereoscopic photogrammetry, interferometric synthetic aperture radar (InSAR), and light detection and ranging (LiDAR), see Wilson [13] for a comprehensive review. Differences between these methods exist in their footprint, cost, resolution and accuracy of the resulting DEM. Whatever method used, however, the resulting DEM will inevitably contain errors that are introduced either during source data acquisition or during data processing. The so-called DEM error hence refers to the difference between the true real world elevations and their DEM representation. Typically, there is a lack of information on the DEM error, which led to notion of 'DEM uncertainty' that refers to what we do not know about the error, see Wechsler [14]. Nowadays, several global DEM databases, e.g. SRTM [15], AW3D30 [16], and TanDEM-X [17], as well as some regional DEM databases [18] are publicly available. Also commercial DEM databases exist that are associated with significant costs, e.g. Hawker et al. [19]. Online initiatives such as OpenTopography facilitate community access and aim at democratizing online availability of high-resolution topography data acquired with LiDAR and other technologies [20]. Despite the broad variety of existing DEM sources, however, we are still facing (and will face in the near future) a very limited availability of high-accuracy DEMs for some regions that are particularly prone to landslide hazards, e.g. in Asia [2]. Whenever using DEM data for simulation-based landslide hazard analysis, it is hence important to be aware of DEM error and uncertainty, and to consider its potential impact on computational run-out analyses and related computational risk assessments. DEM error has arisen researchers' attention since long. Many efforts have for instance been put into quantifying the error associated with specific DEM sources based on data of higher accuracy, e.g. acquired by satellite measurements [21,22], medium footprint LiDAR [23], or GPS survey [15,17,24,25,26]. Meanwhile, a variety of methods have been devised to classify DEM error into various categories [27,28,29]. Due to the complexity of potential influencing factors (sensor technology, retrieval algorithms, data processing, land cover and surface morphology, terrain attributes, see [13,29,30]), these methods can only constrain the DEM error, and will not deterministically correct for it at all grid points. Hence, DEM uncertainty remains, and has to be accounted for in any subsequent analysis that relies on the DEM data. In this circumstance, a stochastic simulation is an effective computational approach to deal with the situation [31]. Instead of considering a single (assumed as accurate) DEM, the fundamental idea of a stochastic simulation in the context of DEM uncertainty propagation is to separate the DEM into a known deterministic contribution and an unknown DEM error. DEM uncertainty is then accounted for by treating the DEM error as a random field consisting of a collection of random variables defined at selected grid points. An ensemble of equiprobable realizations of the random field is then generated based on certain assumptions and available information of DEM error. This could for instance be the so-called root mean square error (RMSE), a minimalistic indicator for the overall error magnitude or a semivariogram that informs about the spatial autocorrelation of the DEM error. Adding the DEM error realizations to the known deterministic DEM contribution results in an ensemble of equiprobable DEM realizations, which can finally be used for a DEM uncertainty propagation analysis. Stochastic simulation methods for DEM uncertainty propagation analyses have been developed since the 1990s and are by now widely applied in many fields, including terrain analysis [31,32,33], flood modeling [19,34,35], soil erosion modeling [36], landslide susceptibility mapping [37], dry block and ash flow modeling [38], etc. With respect to rapid, flow-like landslide run-out modeling, very little work has been done to assess the potential impact of DEM uncertainty, most likely due to the complexity, and hence level of sophistication of the associated process models. Meanwhile, however, advances in computing technology led to computationally feasible and well-developed landslide run-out simulation tools. As one of the most important inputs for these tools, a DEM determines the landslide's flow path. A natural next step is hence to consider the impact of DEM uncertainty in these models, as overlooking DEM uncertainty may lead to a bias of risk management decisions in a wrong direction. The major aim of this study is therefore to describe two different approaches in order to incorporate DEM uncertainty into computational landslide run-out analyses, and to investigate and discuss their feasibility, as well as whether DEM uncertainty is critical to landslide run-out and affects its results. This paper is organized as follows: In section 2, we briefly describe the landslide run-out model used in this study, which is a continuum-mechanical shallow flow model based on the Voellmy-Salm rheology. In section 3, we recall on various methods to account for DEM uncertainty with a major focus on two approaches, namely an unconditional and a conditional stochastic simulation method. The rest of the paper is devoted to investigating DEM uncertainty propagation for rapid, flow-like landslides based on an integrated workflow that combines the aforementioned computational process model (section 2) with the stochastic DEM simulations (section 3). Note, that while in our particular study we chose to use a continuum-mechanical shallow flow process model based on the Voellmy-Salm rheology, the workflow itself is modular and non-intrusive. It would hence also possible to couple the stochastic DEM simulation with any other (DEM based) computational landslide model. Section 4 describes the modular Python-based workflow that we developed in order to set-up and manage the workflow and to interpret its simulation results. We present a series of computational scenarios based upon a historic landslide event in section 5. All scenarios compare the unconditional and conditional stochastic DEM simulation. Finally, section 6 is devoted to a discussion of our results. Important conclusions are drawn in section 7.

Landslide process model
Let {X, Y, Z} denote a fixed Cartesian coordinate system, in which X and Y are the horizontal axes and Z is the vertical axis. The coordinates of a point in the Cartesian coordinate system are denoted by (X, Y, Z) and a topography can then be expressed as the surface mapping of horizontal X and Y coordinates and represents the elevation at each point, namely Z(X, Y ). The mapped topography induces a surface coordinate system {x, y, z}, in which x and y denote tangential directions and z is perpendicular to the surface. As mentioned in the introduction, a variety of numerical landslide run-out models have been developed, among these is the family of depth-integrated shallow flow type landslide models. The latter can be further classified based on their applied basal rheology, e.g. Voellmy, Bingham, Quadratic resistance model, etc., see also [39,40]. Our study relies on a Voellmy-Salm (VS) model, hence a depth-averaged continuum mechanical model using the Voellmy basal rheology. Such models are not only applied to rapid landslides, but also to snow avalanches and certain types of rock fall. In our study the process model will be adopted to the historic landslide event in section 5. The process model along with its underlying assumptions is described in [41,5]. It assesses the slide's dynamics in terms of flow height and depth-averaged velocity, both of which depend on time t and spatial coordinates x and y. They are denoted by H(x, y, t) and U (x, y, t) := (U x (x, y, t), U y (x, y, t)) T re-spectively. The governing partial differential equations are derived from first principles of mass and momentum conservation and read The first equation of system (1) denotes the mass balance, in which H denotes the landslide's height, U x and U y its surface tangential velocity components, andQ(x, y, t) stands for a mass production source term that accounts for erosion of material along the way. Second and third equations denote the x and y momentum balance, in which g x , g y , and g z are the three components of gravitational acceleration vector g, n x and n y are x and y component of the unit vector that opposes the velocity, and µ and ξ stand for dry Coulomb and 'turbulent' friction coefficients respectively. The two friction parameters are mostly determined by back-analysis based on historic events. Finally, c x and c y are velocity shape factors and k a/p denotes the so-called earthpressure coefficient. Both are difficult to determine for realistic slides, yet have been shown to not critically affect the slide's simulation [5]. In this study, we therefore set c x = c y = 1 assuming almost plug flow, and k a/p = 1 since non-hydrostatic internal stresses play a negligible role for flow-like, fully saturated flows [42]. The topographic surface Z(X, Y ) enters the governing equations of the process model implicitly in terms of the spatially varying gravitational acceleration vector g = (g x , g y , g z ) T . Any error and uncertainty present in the topography representation hence also influences the landslide run-out simulation result. The VS model had been first proposed to model snow avalanche [43]. Nowadays, it has been widely applied to other types of gravity-driven rapid mass movements including flow-like landslides [11,44,45,46]. In this study, the proprietary mass flow simulation platform RAMMS [5] which provides a GIS integrated implementation of the VS model is used for landslide run-out modeling. It is integrated as a module of our workflow (see section 4) that is developed for the purpose of DEM uncertainty propagation analysis. For a detailed discussion of the VS model and its implementation, we refer to [5,41].

Methods to represent DEM uncertainty
As before in section 2, we express the topographic surface as a function Z(X, Y ) parametrized in horizontal coordinates X and Y . In practice, the function Z(X, Y ) is often constructed from discrete gridded raster data. We hence assume that a domain of interest D is discretized into the horizontal X and Y direction, which results in a spatial grid defined as The elevation data associated with each grid point D ij is defined as The elevation Z mn of a common DEM data product might be erroneous with respect to the true values as discussed in the introduction. If we denote the true elevation as the DEM error can be expressed as If we knew the error mn , we would be able to recover the real world topographic surface Z * mn . The fact, however, that the error is unknown, or we only have limited information about the error implies an uncertainty to the input of our landslide process simulation. Within this study, we will refer to the uncertainty associated to the unknown DEM error as DEM uncertainty. In this circumstance, each ij is treated as a random variable and mn is accordingly treated as a random field, which consists of a collection of random variables ij . By generating multiple realizations of the random field mn , DEM uncertainty can be represented. This process is widely known as stochastic simulation. It requires a suitable model to describe the jointed uncertainty of all ij based on limited available information of DEM error. The task can be further divided into determining: • the probability distribution function (pdf) of each ij which quantifies local uncertainty at each grid point; • the correlation between different ij which is usually termed as spatial autocorrelation of DEM error.
According to available information on the DEM error, existing approaches that could be used to solve the two issues can be generally classified into two groups: A) unconditional stochastic simulation (USS), and B) conditional stochastic simulation (CSS).
More specifically, USS is only informed with properties of DEM error, e.g. the RMSE, and thus does not honour any actual DEM error values. CSS is informed with certain number of actual DEM error values at reference locations within the DEM, e.g. obtained from higher accurate reference data, and thus could directly honour the actual DEM error values at reference locations [29].

Unconditional stochastic simulation (USS) based on the RMSE
Typically available information about the DEM error provided by DEM vendors is the root mean square error (RMSE). For a set of K reference locations, it is defined as Here, Z * KK = {Z * kk = Z * (X k , Y k ) | (X k , Y k ) ∈ D; k = 1, 2, ..., K} denotes higher accurate elevation values measured at the reference locations and Z KK = {Z kk = Z(X k , Y k ) | (X k , Y k ) ∈ D; k = 1, 2, ..., K} denotes corresponding elevation values based on the DEM. It should be noted that while the RMSE is typically available, this is not true for the reference elevation values Z * KK itself. As stated numerous times in the literature, it is critical that the RMSE only provides a global indication of DEM error magnitude without any information about its spatial autocorrelation. Still, it is by far the most widely used DEM error indicator for many DEM databases and mostly the only available information coming along with DEM products. In this circumstance, USS could be used to represent DEM uncertainty and study its propagation into landslide run-out analysis. Modeling DEM uncertainty based on USS assumes that all local error values ij are independent and fulfill the same univariate Gaussian distribution with a mean (µ) of zero and a standard deviation (σ) equivalent to the given RMSE. Under this assumption, an ensemble of spatially uncorrelated realizations of the random field mn can be generated by randomly assigning error values to each ij according to its local Gaussian probability distribution. In the next step, we have to account for the (unknown) spatial autocorrelation of mn . Potential methods that could be applied are simulated annealing, spatial autoregressive modeling, spatial moving averages, etc., see [14]. Simulated annealing is generally computationally intensive and spatial autoregressive modeling becomes impractical for simulation of large areas [47]. In this study, we use the spatial moving averages method that increases the spatial autocorrelation by filtering spatially uncorrelated realizations with a distance-weighted filter proposed by Wechsler and Kroll [48]. For ij at one grid point of an uncorrelated realization, its value is replaced by the weighted average of ij at all grid points within the filter kernel. The weight decreases with increasing of the distance to the grid point, which is similar to semivariogram trends [48]. The size of the filter denoted as d depends on the maximum autocorrelation length of mn which again is unknown if the RMSE is the only available information. In practice, d is often determined based on the maximum autocorrelation length of the original DEM [14,36]. Though it relies on some assumptions, such as an appropriate choice of correlation length d, the sketched approach is generally applicable if RMSE is the only available information. It may become critical if a DEM contains a systematic bias which means that the mean of Z * kk − Z kk deviates largely from zero. More specifically, if we follow [29,17] in defining mean µ and standard deviation σ as we can express the RMSE in terms of µ and σ as If the number of reference points K is relatively large, (K − 1)/K is close to one. Eq. (8) then indicates that the RMSE is larger than the standard deviation σ if the mean µ deviates from zero. The difference between the RMSE and σ increases with µ increasing. For example, the µ, σ, and RMSE . This means that assuming the standard deviation of the DEM error to be given as the RMSE consequently overestimates the variability of the DEM error if the mean deviates largely from zero. The implications of both issues, namely the fact that the filter size d is unknown and has to be subjectively chosen, and that the RMSE provides an insufficient representation of the DEM error, are investigated in the following study. For convenience, the two issues are referred to as: • unrepresentative RMSE, • subjective d.

Conditional stochastic simulation (CSS) based on higher accurate reference data
This approach requires the availability of higher accurate reference data at certain reference locations, e.g. from higher accurate DEM products, or GPS surveys. Note, that although these data might be subject to error themselves, it is fair to assume this error to be much smaller. This justifies to use the higher accurate reference data as true elevation values Z * KK . Based on Z * KK , we could determine the statistics of the DEM error, e.g. the RMSE, the µ and the σ as discussed in section 3.1. Likewise, we can assess the spatial autocorrelation of the DEM error, e.g. in the form of a semivariogram model (see figure 4). In addition, we know the DEM error at the reference locations exactly, denoted as * KK = { * kk | k = 1, 2, ..., K}. Yet, we still lack knowledge about the DEM error away from the K reference locations, hence the complete random field mn . In that situation, conditional stochastic simulation (CSS) can be used to simulate, i.e. generate realizations of the random field mn . Many geostatistical methods of conditional simulation could be applied, including sequential simulation algorithms, p-field approach, simulated annealing, etc. [49]. In this study, we apply a sequential Gaussian simulation. It is the most attractive technique for stochastic spatial simulation according to Temme et al. [50] and has been widely utilized in DEM uncertainty propagation analysis [31,36].
The sequential Gaussian simulation sequentially samples each local error ij along a random path that consists of all grid points D ij under the multi-Gaussian assumption. This means that assuming the random field mn to satisfy a multivariate Gaussian distribution, hence each ij fulfills a univariate Gaussian distribution denoted as N (µ ij , σ ij ). The essential idea now is that the mean µ ij and standard deviation σ ij are determined sequentially by means of simple kriging based on: the semivariogram model of DEM error that provides covariances in simple kriging equations, and the conditioning information including * KK and previously sampled ij . By making each univariate Gaussian distribution conditional not only to * KK but also to all previously sampled ij , the semivariogram model of DEM error is reproduced in realizations of mn [49]. The process to generate one realization of mn is as follows: 1) determine a semivariogram model to represent the spatial autocorrelation of DEM error based on normal-score-transformed * KK ; 2) define a random path visiting each D ij once; 3) at each D ij , determine N (µ ij , σ ij ) using simple kriging based on the semivariogram model and normal-score-transformed * KK ; 4) sample a value from N (µ ij , σ ij ), assign it to ij , and add ij into normalscore-transformed * KK ; 5) repeat steps 3) and 4) until all D ij along the path are visited; 6) back-transform all sampled ij to the original distribution of * KK . Multiple realizations can be generated by defining different random paths.

Implementation
Studying the impact of DEM uncertainty on landslide run-out modeling is computationally intensive and technically demanding. It includes representing DEM uncertainty in terms of a large number of DEM realizations, conducting numerous landslide run-out modeling based on the DEM realizations, and postprocessing extensive output data. In addition, understanding how DEM uncertainty affects terrain attributes may facilitate us to interpret its impact on landslide run-out modeling. This requires the ability to calculate terrain attributes, e.g. slope, ruggedness, etc. of the original DEM as well as the generated DEM realizations. In this study, we propose a workflow that integrates our own Python implementation of selected aspects of the workflow and existing software as well as toolboxes to solve above mentioned tasks. It is part of our PSI-slide (Predictive Simulation of Slides) package in development that is designed for the purpose of systematically investigating the impact of the various sources of uncertainty on simulating gravity-driven mass movements [51,52]. Herein, we focus on DEM uncertainty. Figure 1 illustrates the workflow. It consists of four modules: 1) DEM uncertainty representation. In this module, we generate an ensemble of N equally probable DEM realizations to represent DEM uncertainty based on available information about DEM error. USS as introduced in section 3.1 is implemented without third party software (USS.py) for cases in which only the RMSE is available. For cases in which higher accurate reference data is provided, CSS as introduced in section 3.2 is implemented by integrating the sequential Gaussian simulation algorithm of the Stanford Geostatistical Modeling Software (SGeMS) [53] into our workflow (SGeMS.py). 2) Landslide run-out modeling. This module is used to conduct N landslide run-out simulations based on the N DEM realizations generated in module 1). In this study we employ the proprietary mass flow simulation platform RAMMS [5] which provides a GIS integrated implementation of the VS model. First, a Python script named SetInput.py is called to set up required inputs for each simulation run. Then a Python script named RAMMS.py starts parallel runs of RAMMS using the Python Scoop module. In the end, a Python script named ExtractOutput.py is called to extract user-specified outputs. 3) Statistical analysis and visualization. This module is used to conduct statistical analysis on the user-specified outputs from module 2) and visualize results. It is mainly based on the Python Numpy and Matplotlib modules. For example, probabilistic hazard map can be produced to indicate potential hazard area. 4) Terrain analysis. This module is used to analyze terrain characteristics of the original DEM and DEM realizations from module 1), which may help us to interpret outputs from module 3). This is achieved by integrating several terrain analysis tools from WhiteboxTools [54] like calculating slope, aspect, ruggedness index, etc. into our workflow (WhiteboxTools.py). It is part of our PSI-slide package in development that is designed for the purpose of systematically investigating the impact of various sources of uncertainty on simulating gravity-driven mass movements [51,52]. The workflow consists of four modules: 1) DEM uncertainty representation; 2) landslide run-out modeling; 3) statistical analysis and visualization; 4) terrain analysis.

Case study
This study is based upon a historic landslide and two DEM sources. For the purpose of DEM uncertainty propagation analysis, we assume one DEM source to be more accurate than the other and then obtain higher accurate reference data from the more accurate DEM source to assess elevation error of the less accurate DEM source. We design a series of computational scenarios based on the higher accurate reference data to study the impact of DEM uncertainty on landslide process simulation for both the case when only the RMSE is available and the case when higher accurate reference data is available. Additional computational scenarios are designed to study the unrepresentative RMSE and subjective d issues as detailed in section 3.1 in the form of a sensitivity analysis.

Scenario background and DEM sources
The historic landslide happened on June 7 2008 on the hillside above Yu Tung Road in Hong Kong due to an intense rainfall event, see figure 2. It was the largest flow-like landslide out of 19 landslides during that event. Around 3400 m 3 material were mobilized and traveled about 600 m until deposit. The landslide event had a severe infrastructural impact, as it led to closure of westbound lanes of Yu Tung Road for more than two months [55]. The Yu Tung Road landslide also served as a benchmark case for predictive landslide run-out analysis at the second Joint Technical Committee on Natural Slopes and Landslides (JTC1) Workshop on Triggering and Propagation of Rapid Flow-like Landslides in Hong Kong 2018 [11]. Two types of DEM data of the Yu Tung Road area had been the basis for this study:  landslide in the GEO report [55].
In this study, we assume the 2m-DEM to be more accurate than the 5 m resolution HK-DTM. Similar to our consideration at the beginning of section 3.2, the 2m-DEM and 5 m resolution HK-DTM correspond to Z * and Z as defined in section 3. A set of higher accurate reference data Z * KK can hence be determined to provide information to represent uncertainty of the 5 m resolution HK-DTM. It should be noted that the 2m-DEM and 5 m resolution HK-DTM were produced in different time periods. After the 2008 Yu Tung Road landslide, debris-resisting barriers and a road had been built in the area within the red circle and blue rectangle in figure 3 (a) respectively. They are reflected in the HK-DTM but not in the 2m-DEM, which leads to large inconsistency between the two DEMs in that area. Therefore, to avoid unrealistically large error of the HK-DTM, data from the 2m-DEM in that area is excluded from higher accurate reference data Z * KK .

DEM realizations
Here, Sph(·) and Exp(·) denote the basic spherical and exponential semivariogram models [49] and h denotes the horizontal distance between any two locations. A comparison between the experimental semivariance values based on * KK {K = 180} and the parametrized semivariogram model given by eq. (9) can be seen in figure 4. Semivariance is a measure of spatial dependence between DEM error values at two different locations. The continuous semivariogram model is fitted to the experimental semivariance values so as to deduce semivariance values for any possible distance h required by simple kriging [49]. The range of the semivariogram model is 180 m. It indicates the maximum autocorrelation length of the HK-DTM error, on which the size of the spatial moving filter d depends (see section 3.1).

DEM uncertainty scenarios
As mentioned in section 3, DEM users are often restricted to DEM error information in the form of a single RMSE value per data product. Rarely, they have higher accurate reference data. In order to account for both situations, two corresponding 'information levels' are considered in the following study.

Number of DEM realizations
The integrity of a stochastic simulation requires a large number of DEM realizations, while more realizations naturally take many computational resources. Thus one has to find a reasonable compromise. Typically, this can be found through a representative convergence study. Since in our study we address the impact of topographic uncertainty on landslide run-out simulation, we analyze the relative change of topographic attributes with an increasing We define an indicator of the relative change similarly as in [32] to investigate the converging behaviour. Taking slope as an example, for a given number n of HK-DTM realizations, we first calculate the standard deviation of slope at each HK-DTM gird point over the n realizations. The calculated standard deviation values at all grid points constitute a grid of standard deviation values. Then we calculate the standard deviation of the grid of standard deviation values, which leads to a single standard deviation value for the given number n. For each n from 1 to 1000, we can correspondingly calculate a standard deviation value. The same procedure is applied to aspect, ruggedness, and elevation. Figure 5 shows plots of normalized standard deviation of the grid of standard deviation values with respect to the number of HK-DTM realizations for the two situations A) and B). It can be seen that for situation A), aspect levels out first, followed by slope, ruggedness, and elevation. Beyond 500 realizations, there is little change of normalized standard deviations. This indicates adding more realizations has little impact on topographic attributes. For situation B), aspect also levels out first while the rest three show less difference. Compared to A), B) converges faster which indicates CSS requires less DEM realizations than USS. Nevertheless, we set N=500 for the remainder of this study both for USS and CSS. We generate 500 HK-DTM realizations for each scenario input set as listed in table 1.

Statistics of DEM error realizations
In order to conduct a further quality check of our implementation of both USS and CSS, we investigate the corresponding DEM error realizations of the USS N=500 {RMSE=3.3, d=180} and CSS N=500 scenarios, denoted as USS Error  This also corresponds to the assumption underlying CSS that each ij fulfill a univariate Gaussian distribution with a mean µ ij and standard deviation σ ij given by the simple kriging estimate and simple kriging standard deviation at D ij (see section 3.2).

Landslide process simulation setup
With the DEM realizations generated in section 5.2, we can study the impact of DEM uncertainty on landslide process simulation. Here, we introduce the key inputs and our setup for the process simulation.

Model input
Release zone area and fracture height, friction parameters, and a DEM are three key inputs for performing a deterministic landslide process simulation based on the VS model and utilizing the mass movement simulation platform RAMMS [5]. For all scenarios, we consistently use the release zone area as provided for the benchmark exercise during the second JTC1 workshop, which match that of the 2008 Yu Tung Road landslide [11] as shown in figure 7 (b). The fracture height is assumed to be 1.2 m leading to a release volume of around 2900 m 3 based on the 5 m resolution HK-DTM. The friction parameters µ and ξ used in this study are 0.105 and 300 m/s 2 respectively. They are suggested in the GEO report issued by Geotechnical Engineering Office of Hong Kong, which are obtained using back-analysis with information from a video capturing the lower portion of the landslide and detailed field mapping after the landslide [55]. The HK-DTM and all HK-DTM realizations generated in section 5.2 are used as DEM inputs. Entrainment is not considered in this study.

Simulation ensembles
We denote a deterministic landslide process simulation based on a DEM as a simulation run and N deterministic landslide process simulations based on N DEM realizations as a simulation ensemble. The following deterministic simulation and simulation ensembles are conducted based on the original HK-DTM and the aforementioned computational scenarios, see

Results and discussions
This section is organized according to the simulation ensembles introduced in section 5.3.2. Section 6.1 presents the results of the deterministic simulation HK-DTM which serves as the basis for all following discussions. Section 6.2 is devoted to analyze the impact of DEM uncertainty on landslide process simulation in the case of RMSE only (USS N=500 {RMSE=3.3, d=180} ensemble) and available higher accurate reference data (CSS N=500 ensemble). In section 6.3, the unrepresentative RMSE and subjective d issues are discussed based on the ensembles described in section 5.3.2 4) and 5).

Deterministic simulation HK-DTM
In a continuum mechanical landslide process model such as used for this study and introduced in section 2, the landslide flow behaviour is characterised by its spatially varying height and velocity distribution over time, denoted as H(x, y, t) and U (x, y, t). For the purpose of landslide hazard assessment and mitigation measure development, hence maximum height and velocity data through the duration of the landslide are most informative. Thus, we focus on the maximum values of H(x, y, t) and U (x, y, t) over all time, denoted as Figure 7 (a) and (b) show H max (x, y) and U max (x, y) as given by the deterministic simulation HK-DTM. It should be noted that there is a relatively high elevation area at the end part of the channel in the HK-DTM as denoted within the red circle in figure 7 (b). It corresponds to the construction of debris-resisting barriers after the 2008 Yu Tung Road landslide as introduced in section 5.1. The flow material is decelerated and held back here. We will come back to this point latter in section 6.2.1. Landslide run-out distance is often characterised in terms of its apparent friction angle. The tangent of the apparent friction angle is equal to the ratio of the landslide fall height and the run-out distance [56]. The apparent friction angle evaluated from the deterministic simulation is 16.80 • . These results are used as reference to assess the impact of DEM uncertainty in the following discussions.

USS N=500 {RMSE=3.3, d=180} ensemble and CSS N=500 ensemble
While it is straightforward to present results of a deterministic simulation run as shown in section 6.1, a stochastic simulation based ensembles of N simulation run call for tailored statistic to manage and interpret the extensive output data. First, we define the hazard probability P (x l ,y l ) at one location (x l , y l ) as the frequency of H max (x l , y l ) exceeding a certain pre-defined height threshold value, hence where H n max (x l , y l ) denotes the maximum flow height at location (x l , y l ) for the n-th simulation run of the corresponding ensemble. P n (x l ,y l ) hence informs whether location (x l , y l ) is within the hazard area of the n-th simulation run for a given threshold, and P (x l ,y l ) about the resulting hazard probability at location (x l , y l ) considering the complete ensemble. Here, the threshold is set as 0.1 m which matches the cut-off threshold of the deterministic simulation HK-DTM in figure 7 (a). Evaluation of hazard probabilities at all locations then gives rise to a probabilistic hazard map [38], which provides an overall view of the DEM uncertainty impact. Besides assessing the overall impact of DEM uncertainty in terms of the probabilistic hazard map, we will also discuss the local impact of DEM uncertainty on dynamic flow properties, focusing on H max (x, y) and U max (x, y) at locations along the channel bottom and the channel cross section denoted in figure 7 (b). Figure 8 (a) and (c) show the probabilistic hazard map for both USS N=500 {RMSE=3.3, d=180} ensemble and CSS N=500 ensemble. It can be seen that the potential hazard area is much larger than the deterministic hazard area for both ensembles. The difference between the deterministic and the ensemble based hazard area is most pronounced in area 1-3 for USS N=500 {RMSE= 3.3, d=180} ensemble and in area 3 for CSS N=500 ensemble. Figure 8  compared to other parts of channel banks. Flow material could be more easily diverted to area 1 and area 2 where elevations are relatively low and some local 'small channels' exist. Area 3 could also be regarded as less 'well defined' since it is relatively flat and thus is sensitive to DEM uncertainty [50].

Probabilistic hazard maps
The change of each topographic characteristic has corresponding impact on landslide run-out behaviour. Specifically, • if banks of the deterministic channel were dampened out in DEM realizations, flow material tends to spread out along channel cross section direction and travel distance would be shorter; • if the relatively high elevation area that holds back flow material was dampened out, flow material tends to travel further; • increasing topographic roughness leads to higher simulated momentum losses and shorter travel distance as pointed out by McDougall [10]. The black outline plotted on the hazard maps represents the deterministic hazard area (see figure 7 (a)). In the boxplots, the blue star denotes the apparent friction angle of the deterministic simulation HK-DTM (see section 6.1). The difference between the deterministic and the ensemble based hazard area is most pronounced in area 1-3 for USS N=500 {RMSE=3.3, d=180} ensemble and in area 3 for CSS N=500 ensemble. Our main findings are: 1) accounting for DEM uncertainty may significantly increase the potential hazard area; 2) the potential hazard area is highly related to the variability of DEM error and topographic details of the original DEM; 3) USS based on the RMSE only may overestimate the spread of potential hazard area and underestimate the travel distance due to an unrepresentative RMSE (e.g. not bias-corrected) that overestimates the variability of DEM error.
For USS N=500 {RMSE=3.3, d=180} ensemble, the variability of DEM error is relatively large, e.g. around 3.3 m governed by the not-bias-corrected RMSE based on * KK {K = 180} (see figure 6 (c)). In this situation, both the north bank near area 1 and south bank near area 2 as well as the relatively high elevation area at the end part of the channel are possible to be dampened out in HK-DTM realizations. For CSS N=500 ensemble, the variability of DEM error is relatively small, e.g. around 1.5 m governed by the standard deviation (σ) based on * KK {K = 180} (see figure 6 (d)). In this situation, the banks tend to remain 'well defined' while the relatively high elevation area is possible to be dampened out in HK-DTM realizations. Thus, area 1 and area 2 are possibly subject to hazard in USS N=500 {RMSE=3.3, d=180} ensemble but less likely in CSS N=500 ensemble. As mentioned above, area 3 is a flat area which is sensitive to DEM uncertainty. Furthermore, it locates near the deposition, around which the impact of upstream DEM uncertainty seems to accumulate. Thus, it is highly affected in both ensembles. The apparent friction angle distribution is determined by a combined effect of change of channel banks, change of the relatively high elevation area at the end part of the channel, and increasing topographic roughness. For USS N=500 {RMSE=3.3, d=180} ensemble, deteriorated channel bank representation and increasing topographic roughness make flow material travel shorter distance, e.g. larger apparent friction angle, while deteriorated relatively high elevation area representation allows flow material to travel further, e.g. smaller apparent friction angle. For CSS N=500 ensemble, channel banks are likely to remain 'well defined' and the degree of topographic roughness increase is lower due to its relatively small variability of DEM error compared to USS N=500 {RMSE=3.3, d=180} ensemble. Thus, flow material in CSS N=500 ensemble tends to travel longer distance, e.g. smaller apparent friction angle, compared to USS N=500 {RMSE=3.3, d=180} ensemble. In summary, we can conclude from the probabilistic hazard maps and boxplots of apparent friction angle distribution that: • accounting for DEM uncertainty may significantly increase the potential hazard area; • the potential hazard area is highly related to the variability of DEM error and topographic characteristics of the original DEM; • USS based on the RMSE only may overestimate the spread of potential hazard area and underestimate travel distance due to a not-biascorrected RMSE that overestimates the variability of DEM error.
It should be noted that the probabilistic hazard map here is constructed based on maximum height and a pre-defined threshold. In simulation-based hazard assessment practice, one may indicate potential hazard using other indicators, e.g. maximum momentum that reflects the impact pressure, etc. and correspondingly modify the threshold value. In this case, our workflow is easily extendible to account for other indicators.

Dynamic flow properties
The left column in figure 9 shows elevation, maximum height and maximum velocity at locations along the channel bottom based on USS N=500 {RMSE=3.3, d=180} ensemble. It is evident that both maximum height and maximum velocity at these locations largely vary from that of the deterministic simulation. Specifically, the mean of maximum height (maximum velocity) values at all the locations based on the deterministic simulation is 1.28 m   The above observations result from similar factors as discussed in section 6.2.1. Due to DEM uncertainty, • ensemble-based flow dynamic properties are likely to vary from that of the deterministic simulation. Larger variability of DEM error is likely to result in more extreme values. As discussed in section 6.2.1, the variability of DEM error for USS N=500 {RMSE=3.3, d=180} ensemble is larger than that for CSS N=500 ensemble due to unrepresentative RMSE issue. Thus the variation range of USS N=500 {RMSE=3.3, d=180} ensemble-based flow dynamic properties is generally larger than that of CSS N=500 ensemble-based flow dynamic properties, e.g. larger mean of ensemble-based 90% confidence interval (the trend would be more clear if we also consider outliers outside 90% confidence interval). • banks of the deterministic channel may be dampened out in DEM realizations. Deteriorated channel bank representation makes flow material more spread out along channel cross section direction. This could lead to smaller ensemble-based mean of flow dynamic properties at channel bottom locations, compared to flow dynamic properties of the deterministic simulation. It can be directly seen in figure 10, which displays results of one channel cross section. Also, due to larger variability of DEM error, flow material in USS N=500 {RMSE=3.3, d=180} ensemble is more spread along channel cross section direction, resulting in smaller ensemble-based mean of flow dynamic properties at channel bottom locations compared to CSS N=500 ensemble. This can also be seen in figure 10. • topographic roughness in DEM realizations tends to increase. As pointed out in section 6.2.1, increasing topographic roughness results in higher simulated momentum losses and thus smaller flow dynamic properties on average. The higher the degree of topographic roughness increase, the higher the simulated momentum losses and the smaller the flow dynamic properties. This also contributes to smaller ensemble-based mean of flow dynamic properties at channel bottom locations, compared to flow dynamic properties of the deterministic simulation, as well as smaller USS N=500 {RMSE=3.3, d=180} ensemble-based mean of flow dynamic properties at channel bottom locations, compared to CSS N=500 ensemble.
Based on the ensembles' dynamic flow properties we can conclude that: • accounting for DEM uncertainty may significantly affect dynamic flow properties, e.g. maximum height and maximum velocity, hence any hazard assessment that is based on landslide dynamics; • USS based on the RMSE only may overestimate the range of dynamic flow properties and underestimate ensemble-based mean of dynamic flow properties due to an unrepresentative RMSE that overestimates the variability of DEM error. ensembles. This indicates that whether spatial autocorrelation is considered or not may make a difference but the extent of spatial autocorrelation has less influence on simulation results. As we know spatial autocorrelation to be present in topographic data but often lack information on its exact autocorrelation length, this is actually good news for practical hazard assessment studies. Comparing subfigures 11-14 (a) with (c), it can furthermore be seen that the results of the USS N=500 {RMSE=1.5, d=180} ensemble are quite close to the results of the CSS N=500 ensemble. The USS N=500 {RMSE=1.5, d=180} ensemble is informed with the bias-corrected RMSE (namely the true standard deviation, in our case 1.5 m, see figure 3 (b)). It indicates that if a bias-corrected RMSE is given, USS is possible to provide reasonable results considering the extent of spatial autocorrelation has less influence on simulation results.    All in all, we find that: • the results of USS are in general more sensitive to values of the RMSE and less sensitive to values of d; • an unrepresentative RMSE that overestimates the variability of DEM error may overestimate the potential hazard area and extreme values of dynamic flow properties; • whether or not spatial autocorrelation of DEM error is considered can make a difference of ensemble-based simulation results; • if a bias-corrected RMSE is given, it is possible to obtain reasonable ensemble-based results using USS.

Conclusions
In this paper, we investigated different approaches to study the impact of topographic uncertainty on simulation-based flow-like landslide run-out analyses. Based upon a historic landslide event for which two DEM data sets of different accuracy had been available, we presented a series of computational scenarios. Unconditional and conditional stochastic simulation are conducted to generate DEM realizations, both for the case in which only the RMSE is available, and for the case in which reference data of higher accuracy is available. The computational workflow including the stochastic simulation to generate the DEM realizations and the landslide run-out simulation is implemented as a modular Python-based package. How topographic uncertainty propagates into results of landslide run-out analysis is discussed in detail. In addition, we addressed the two major issues of purely RMSE-based unconditional stochastic simulation, e.g. the fact that not-bias-corrected RMSE overestimates the variability of DEM error (referred to as unrepresentative RMSE in our study) and the fact that determining the size of the spatially moving filter in the absence of further information on the spatial DEM error structure is often subjective (referred to as subjective d in our study). Our main conclusions are: • DEM uncertainty significantly affects simulation-based landslide runout modeling depending on how well the underlying flow path is represented, which is determined by topographic characteristics of the original DEM and the variability of DEM error. For the same degree of variability of DEM error, the less 'well defined' parts of the flow path in the original DEM are more likely to be affected and thus leads to change of flow behaviour at these parts. Also, an increasing variability of DEM error leads to an increased impact. More specifically, with increasing variability of the DEM error, the potential hazard area and extreme values of dynamic flow properties are likely to increase. This shows the importance of considering topography induced uncertainty for simulation-based landslide hazard assessment rather than simply trusting results of a deterministic simulation if a high accuracy of DEM source is not guaranteed. Also, a preliminary terrain analysis may give some hints on areas that will potential be affected by a topographic uncertainty study. • Both unconditional and conditional stochastic simulation methods can be applied to study DEM uncertainty propagation in landslide run-out modeling. Their main difference is that the computationally performant unconditional stochastic simulation can be conducted based on RMSE information only, while the computationally costly conditional stochastic simulation requires the availability of higher accurate reference data and is thus more accurate. The higher accurate reference data provides additional information on the DEM error structure, e.g. its spatial autocorrelation. If the DEM does not contain systematic bias or alternatively a bias-corrected RMSE is provided, the unconditional stochastic simulation yields reasonable results. Otherwise, the assumptions underlying the unconditional stochastic simulation lead to an overestimation of the DEM error variability, which leads to an increased potential impact of DEM uncertainty on the potential hazard area and extreme values of dynamic flow properties. In particular, our study shows that if no higher accurate reference data is available or if computational costs of an conditional stochastic simulation would be too large, the results of a RMSE-based unconditional stochastic simulation can still be used to provide an upper bound of the potential hazard area as well as extreme values of flow dynamic properties for a hazard assessment to take topographic uncertainties into account. • Results of an unconditional stochastic simulation are in general sensitive to the RMSE value as well as sensitive to the fact whether or not the DEM error's spatial autocorrelation is considered. If the latter is taken into account, results are less sensitive to actual value of the DEM error's maximum autocorrelation length. This indicates that determining a representative RMSE may be more important than finding a correct maximum autocorrelation length, while the DEM error's spatial autocorrelation should not be ignored for simulation-based landslide hazard assessment.