## Abstract

This paper describes a new methodology to quantify the variation in the output of a computational fluid dynamics model for block and ash flows, when the digital elevation model (DEM) of the terrain and other inputs are given as a range of possible values with a prescribed uncertainty. Integrating these variations in the possible flows as a function of input uncertainties provides well-defined hazard probabilities at specific locations, i.e. a hazard map. Earlier work provided a methodology for assessing hazards based on variations in flow initiation and friction parameters. This paper extends this approach to include the effect of terrain error and uncertainty. The results are based on potential flows at Mammoth Mountain, CA, and Galeras Volcano, Colombia. The analysis establishes the soundness of the approach and the effect of including the uncertainty in DEMs in the construction of probabilistic hazard maps.

## 1. Introduction

Perhaps, the most fundamental product created by field volcanologists to characterize the potential for destruction of a volcano is the hazards map. Often a reasonable hazards map can be made when the distribution of deposits of a given type are well exposed, and easily dated and mapped. In general, however, difficult logistics or a paucity of previous work may render understanding of a volcano's history quite incomplete. Moreover, the depositional record on the flanks of a volcano cannot often be assumed to be very complete.

Several studies have explored the use of computational fluid dynamics (CFD) models to produce volcanic hazard maps for a variety of phenomena at a number of volcanoes (Hooper *et al.* 2003; Stinton *et al.* 2006; Murcia *et al.* 2010; Procter *et al.* 2010; Sheridan *et al.* 2010). Hazard maps for ground-hugging flows that are constrained by the terrain, such as pyroclastic density currents and lava flows, are often constructed using a digital representation of the terrain (Takahashi & Tsujimoto 2000; Dalbey *et al.* 2008). Usually, these terrain representations are digital elevation models (DEMs). For this type of study, terrain elevation is rightly recognized as the most essential and fundamental of variables in geographical analysis (Mitasova *et al.* 1996; Atkinson 2002; Wechsler & Kroll 2006; Stefanescu *et al.*, in press). Dalbey *et al.* (2008) introduced procedures for constructing hazard maps using ensembles of CFD model simulations (the TITAN2D code; Patra *et al.* (2005)) of such flows constructed by establishing probability distributions of input uncertainties in flow initiation (location and volumes) and by sampling them. The important contribution of DEM uncertainty to the variability of the flow outcomes was not included in that work as there were no procedures readily available. This work is focused on addressing this lacuna.

A digital representation of a terrain surface is an approximation of reality and is often subject to significant error (Mitasova *et al.* 1996). The error is usually not known in terms of both magnitude and spatial distribution. There are in fact large uncertainties associated with the construction of DEMs. Wechsler & Kroll (2006) showed that DEMs contain errors derived from a variety of sources, such as sampling, measurement and interpolation, and these errors cannot always be well estimated. When such DEMs with errors are used in *a posteriori* analysis, such as in simulations of flows, the errors propagate to the predicted flow.

The most important part of DEM error propagation analysis is the appropriate characterization of the error within the DEM itself, including information about its distribution and spatial structure (Shortridge 2001). DEM vendors generally provide users with a measure of vertical accuracy in the form of the root mean squared error (r.m.s.e.) statistic. However, many papers have reported on the limitations of a single value of accuracy, stressing that DEM error is spatially variable and highly correlated (Wechsler & Kroll 2006; Darnell *et al.* 2008). Also the magnitude of the DEM error is closely related to the characteristics of the terrain surface. For example, slope will influence interpolation procedures.

DEM error propagation analysis was introduced to the geographic information system (GIS) community in the early 1990s. In the work of Heuvelink *et al.* (1990), error propagation in calculating slope and aspect was represented using Monte Carlo simulation. It was shown that standard deviations of slope and aspect were higher than expected. The effect of error in the DEMs on the erosion models was emphasized. A method used by Weng (2002) in quantification of the uncertainty of DEMs was to create various DEMs using different interpolation methods and to examine the r.m.s.e. from the source map, sampling and measurement error, and the interpolation process. It was concluded that the r.m.s.e. can be used as a general indicator of DEM uncertainty. In the literature, DEM error without spatial autocorrelation was considered to be a worst case scenario (Heuvelink *et al.* 1989; Van Niel *et al.* 2004; Oksanen 2006), but no analysis based on terrain morphology and the effect of different DEMs was done. Wechsler & Kroll (2006) developed four different methods for representing the spatial dependence of error through random fields to assess the effect on topographic parameters of the DEM uncertainty. The study showed that uncertainty in the DEM is manifested at higher elevations in locally steeper slopes, on both slope and elevation maps. Florinsky (1998) showed that the effect of DEM uncertainty on the accuracy of slope and aspect estimation cannot be determined by using data from topographic maps or field surveys, because accurate derivatives cannot be determined.

One key feature of spatial data is the autocorrelation of observations in space. Generally, spatial autocorrelation refers to the correlation between the same attribute at two locations. Observations in close spatial proximity tend to be more related than are observations at larger distance or separation. Errors in spatial data (such as incorrect elevation values assigned to a point) are spatially autocorrelated. The effect of correlated DEM error has been investigated in the literature (Fisher 1991; Goodchild *et al.* 1992). It was shown that not only is error spatially variable throughout a DEM, but, within the elevation model, the error value of an individual grid cell is related to the error in neighbouring cells. Unfortunately, DEM providers do not include information regarding the spatial dependence or spatial relationship of errors.

Stochastic modelling uses stochastic conditional simulation to generate multiple equally likely representations of an actual terrain surface. Ehlschlaeger & Shortridge (1996) and Hunter & Goodchild (1997) computed a normal distribution of maps or realizations to reproduce the spatial autocorrelation encountered in the original error surface, filtered using a Gaussian convolution filter, with kernel sizes derived from autocorrelation analysis of the original error surfaces.

Various researchers have applied stochastic techniques to evaluate uncertainty in DEM data. Ehlschlaeger & Shortridge (1996) stochastically simulated error in a DEM to evaluate the impact of DEM uncertainty on a least-cost-path application. Hunter & Goodchild (1997) investigated the effect of simulated changes in elevation at different levels of spatial autocorrelation on slope and aspect calculations. Hebeler & Purves (2008) produced uncertainty surfaces to show the impact of DEM uncertainty on an ice sheet model. Darnell *et al.* (2008) developed a fuzzy framework to examine the probable and possible uncertainties in classifying landslide hazard.

The aim of this paper is to quantify the variation in the output of a computational flow model for block and ash flows, when the model inputs, including the elevation values represented in the DEM, are uncertain or given as a range of possible values. Integrating these variations in the possible flows as a function of input uncertainties provides well-defined data on the probability of hazard at specific locations, i.e. a hazard map (Dalbey *et al.* 2008). In particular, the focus here is on assessing the influence of DEM uncertainties (along with uncertainties in initial size and location of the avalanche, and the internal and bed friction angles). There is uncertainty in all of these inputs, which can be represented using either field data or stochastic methods. The distribution or the range of the parameters can be obtained from laboratory and field instruments for friction angles, and historical records of flow frequency and magnitude for size of the initial failure. Stochastic methods are used to assess the uncertainties in the DEMs: the first method consists of a perturbation of the elevation based on the measured error model, while the second method represents an unconditional stochastic simulation (Ehlschlaeger & Shortridge 1996). Both methods generate multiple likely representations of the actual terrain, while the second one accounts for the spatial autocorrelation between elevation points. The effect of DEM uncertainty and its impact on the model output is analysed by constructing a hazard map and performing a ‘probability analysis’ for two volcanoes with different morphology: Galeras Volcano, Colombia, and Mammoth Mountain, CA. The second approach adapted here is based largely on the method of Ehlschlaeger & Shortridge (1996), which uses the difference between two independent DEMs to train a Gaussian model of error.

The next section presents the basic methodology for generating ensembles of DEMs representative of the true DEM. Subsequent sections summarize the TITAN2D flow simulation tool and its use in a systematic hazard analysis. The hazard analysis tool itself uses ensembles of TITAN2D simulations to construct statistical surrogate models of flow outcomes at different locations as a function of model inputs, such as flow volume and resistance to flow as modelled by a Coulomb frictional law. Sampling of these surrogates leads to the construction of effective hazard maps that reflect the range of uncertainty in the model inputs.

## 2. Methodology

In previous work (Stefanescu *et al.*, in press), the effect of DEM variability on the output of TITAN2D was investigated by comparing an output variable—maximum flow depth over the entire time of simulation—from different DEMs of the same site. These DEMs were obtained from different techniques at different resolution. Two types of analysis were performed: a qualitative analysis and a statistical analysis. The qualitative analysis consisted of a comparison of the footprint of the flow, extended to a pixel-based classification. The pixels were classified into inundated and non-inundated classes. For the statistical analysis, a Kolmogrov–Smirnov test was performed to assess whether the two output datasets differed significantly. The conclusion was that, for moderate and small-scale flows, use of different DEMs affects computation of accurate footprints of the flow.

This conclusion motivated the present study to examine the effect of DEM uncertainty by creating a model of the error and sampling it to create an ensemble of possible terrains. The flow simulation is then run on every member of this ensemble.

Naive, cell-by-cell approaches to treating DEM uncertainty quickly lead to the use of thousands if not millions of random variables, resulting in a computationally infeasible problem. On the other hand, the error model described above can be parametrized with one or two random variables. The parametrization methods are based on the assumption that the available DEM is a representation of the terrain to which errors have been added because of instrumental uncertainty. Therefore, the DEM can be assumed to be one of an infinite number of elevation realizations.

### (a) Method 1

In this paper, two ‘types’ of DEMs are available of each mountain, which are used in creating DEM-to-DEM difference maps. Different realizations of the terrain were constructed by adding to one DEM—considered to represent the ‘true’ elevation—a ‘random’ perturbation. Since any two types of DEMs are obtained using different techniques, the difference between them can be added to that which is assumed to be the ‘true’ DEM to give a set of possible DEMs. Thus, the resulting realizations are consistent with the available set of DEMs. Randomness in the perturbations is created by multiplying the difference map with a scalar random variable *ξ*, which is normally distributed between 0 and 1,
2.1where *R* is a realization of the terrain, *M* is the DEM that best represents the terrain (the ‘true’ DEM) and Diff is the difference map. In this way, we can define a set of DEM realizations using only one random variable.

### (b) Method 2

For elevation, data at any grid point in a DEM tends to be related to data from the nearby points. This is the principal motivation of method 2, based on the work of Ehlschlaeger & Goodchild (1994). Difference maps can be constructed if more than one DEM exists for the same location. Such maps are termed error maps and are generated by subtracting the lower quality DEM from the higher quality DEM (i.e. the ‘true’ DEM). These maps are spatially autocorrelated. Random fields can be used to represent these spatially autocorrelated data points. Let be a continuous random field used to characterize unknown elevation errors (differences). The random field function is implemented in the function *r.random.surface* (Ehlschlaeger & Goodchild 1994) of the Geographic Resources Analysis Support System (GRASS) GIS (Mitasova *et al.* 1996), and generates fields obtained using a normal distribution (mean of 0.0 and variance of 1.0). The random field function derives its spatial dependence from the use of a distance-based decay filter function. The following equation is used to generate the random field:
2.2and
2.3where is the set of potentially influencing points in a given area, , *w*_{u,v} is the spatial autocorrelative effect between points and , *ϵ*_{v} is a Gaussian random variable with a mean of 0 and variance of 1, *d*_{u,v} is the distance between *u* and *v*, *D* is the minimum distance of spatial independence, *E* is the distance decay exponent and *F* is the distance at which errors are completely correlated.

A set of random fields is calibrated to the spatial variation of the field being simulated using a correlogram function. This is done by fitting the correlogram and choosing the best descriptive parameters of the random field (the minimum distance of spatial independence, the correlated distance decay exponent and the filter parameter) in a weighted least-squares estimator implemented in GRASS's *r.lags.difference*. After running hundreds of tests with multiple combinations of *D*, *E* and *F*, the best random field was found by fitting the error map characteristics, such that the sum of least-squares difference between an error field's correlogram and the target correlogram is minimized. Figure 1 shows a sample error map correlogram and several trial correlograms closely fitting it. From equation (2.3), it can be seen that the parameters *D*, *E* and *F* influence the shape/look of the correlogram. It was noted that the main impact of the exponent value is to characterize the roughness of the texture of the random surface. Surface roughness decreases as the exponent value gets closer to 1.0. Once the parameters are set to a certain value as determined above, one is able to sample from a normal distribution value for *ϵ*_{v} as given in equation (2.2) to generate a possible perturbation of the provided DEMs. In this way, a normal distribution of possible terrain maps is produced, where the mean of the distribution represents the original DEM used as the ‘true’ surface.

The correlogram model was used with sequential Gaussian simulation to generate a set of error map realizations. Each error realization was added to the ‘true’ DEM indicated as , to generate equally probable realizations of the topography for the error structure of a DEM under consideration,
2.4where is a realization of the elevation dataset , is a group of sets of spatially uncorrelated sample points in and *ϵ* is a Gaussian random variable with mean 0.0 and variance 1.0. and variance are mean and variance, respectively, of all sets in . specifies the random field as defined in equation (2.2). Hence, this methodology parametrizes the DEM using only two Gaussian random variables, *ϵ*_{v} and *ϵ*.

### (c) Digital elevation model realizations

Many DEM users are aware that DEM uncertainty affects the results of their application; however, in most cases, the DEM is accepted as the true representation of the Earth's surface. In this section, two methods for generating multiple realizations of the terrain are presented for both Galeras Volcano and Mammoth Mountain, to test whether it is safe to assume that the representation of topography is acceptable as it is.

The motivation for creating a process to generate realizations of the DEM was to incorporate the DEM as one of a host of uncertain input parameters for TITAN2D simulations and consequent hazard map calculations. One working hypothesis is that the DEM contributes a significant proportion of the variance to simulated flow, and hence the hazard map output. For sampling the input parameter space, a Latin hypercube sampling (LHS) was implemented. LHS is commonly used in computer sampling experiments (McKay *et al.* 1979; Sacks *et al.* 1989) mainly because it is computationally cheap to generate and can cope with many input variables. This sampling can also have relative small variance when measuring output variance.

For Galeras Volcano, two test DEMs at 30 m spacing were considered for the analysis. The Shuttle Radar Topography Mission (SRTM) 30 m DEM was derived by spline interpolation from a 90 m DEM of southern Colombia using radar data collected in 2000, while the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEM was calculated at the Jet Propulsion Laboratory (JPL) using orthorectified imagery from 12 January 2010 (figure 2*a*). The ASTER dataset was used as a surrogate for the ‘true’ elevation while the SRTM dataset was used in creating the error model.

Two 30 m resolution DEMs derived from independent techniques were used for Mammoth Mountain. A topographic synthetic aperture radar (TOPSAR) dataset was considered to be the ‘true’ elevation, while an SRTM dataset was used in creating the error map. A rectangular area of approximately 42 km^{2} was defined within the TOPSAR and SRTM DEMs (figure 2*b*).

For method 1, 64 DEM realizations were created and used as input parameters for the TITAN2D simulator along with uncertain parameters presented in §3*c*. The input space is defined by seven parameters.

As described above for method 2, realizations of the terrain surface were created by taking into consideration the spatial autocorrelation of the error. The error map was obtained by subtracting the elevation of a given DEM from the ‘true’ elevation at each location. The correlogram for the difference map was calculated to determine the range of spatial dependence of elevation points. It was found that spatial dependence persisted above a threshold value of the correlogram cross-correlation coefficient of 0.4 to a distance of 2.5 km for Galeras and 2.1 km for Mammoth. To determine the probability distribution function for the stochastic simulation, 91 sets of spot locations were selected from the map, each set containing 91 points; all pairs of points were separated by more than 2.5 km or 2.1 km, respectively. For each DEM, probability distribution function statistics were derived. The random field parameters were chosen after testing more than 400 random field parameters for the smallest difference between the error model correlogram and the random field. This occurs when the minimum distance of spatial independence *D*=2500, the distance decay *E*=0.8 and the filter parameter *F*=400 for Galeras and *D*=2100, *E*=0.7 and *F*=350 for Mammoth. A total of 64 equally probable potential elevation surfaces of the area having a 30 m resolution were generated.

### (d) Hazard map construction

There are numerous ways to create a volcanic hazard map based on CFD modelling. The traditional Monte Carlo method can be used if it is assumed that uncertainty in model input parameters is the main restriction to the knowledge of future events at a given volcano. This is the case, for example, if it is known that block and ash flows are common at a given volcano, but it is difficult to know the size or volume of potential future events. Although Monte Carlo is relatively simple to implement, it converges slowly and is unaffordable computationally because of the number of time-consuming simulations. A single TITAN2D run might take 20 min on a single processor. To obtain three-digit accuracy in the expected value of a specified function would require a million runs. One million runs of 20 min calculations running non-stop on a 64 processor would take 217 days (Dalbey *et al.* 2008).

Here, a brief description of the use of a hierarchical emulator that significantly reduces computational cost is presented; a detailed discussion of the methodology can be found in Dalbey (2009) and Dalbey *et al.* (2008). An emulator can be thought of as a fast statistical surrogate for a single numerical model simulation (a simulator). The process of computing a hazard map for block and ash flows with uncertain model inputs introduced by Dalbey (2009) is described. Two-level construction of a group or ensemble of emulators is used to include a separation of uncertain inputs and geographical coordinates. The process starts by identifying the model inputs whose uncertainties will drive the process. In this case, the uncertain flow inputs used are volume and shape, starting location, basal and internal friction angles, and finally topography, as given by the DEM. For the resulting eight-dimensional parameter input space, an LHS was performed to determine parameter values at which simulations were to be run. As priors for the emulator, simulation outputs for each of these input parameter vectors were stored at 64 grid points. The sample size is consistent with other numerical experiments of this type existing in the literature (McKay *et al.* 1979; Sacks *et al.* 1989; Mitasova *et al.* 1996).

The output variable of interest for application in this paper is the field of maximum flow depth over time for each spatial position, at each of the downsampled input parameter grid points. Tesselations of the geographical coordinate space and the parameter input space are constructed (a Delaunay triangulation was used). At a designated location, **x***, of the input parameter plus spatial coordinate space at which the hazard is to be computed, the covering simplex *S**_{x} of the parameter space is identified, and all nodes of that simplex are enumerated, as are all nodes within a neighbourhood (two hops in the tesselation) of the covering simplex nodes. For each such two-hop node, the tesselation was performed in the spatial coordinates followed by an evaluation of all emulators constructed over these nodes. These coordinate space emulators to (the coordinate components of) **x*** by barycentric weighting were averaged; notice that there will be an emulator for each parameter input sample point. Now in the input parameter space, construct a tessellation of the two-hop nodes and average the emulators to **x*** by barycentric weighting of the fine-scale emulator. The emulator is now readily and quickly evaluated for each evaluation. The hazard map construction can now proceed by treating the emulator as a surrogate for the simulator in the classical Monte Carlo procedure. For any point in the domain, it can now be exercised like the simulator to get potential flows and hence exceedance probabilities.

### (e) TITAN2D and flow simulations

TITAN2D (Patra *et al.* 2005; Sheridan *et al.* 2005) was developed for modelling dry geophysical granular flows, such as debris avalanches and block and ash flows. Given a digital elevation map specifying the topography of a volcano and the values of input parameters, including the initial volume of erupted material and the friction angles, TITAN2D calculates the flow depth and velocity at any location throughout the duration of an event. The TITAN2D code combines numerical simulations of a natural granular flow with digital terrain data. It is based on a depth-averaged model for an incompressible granular material governed by Coulomb-type friction interactions (Savage & Hutter 1989). The governing equations are obtained by applying conservation laws to the incompressible continuum, providing appropriate constitutive modelling assumptions, and then taking advantage of the shallowness of the flows (flows are much longer and wider than they are deep) to obtain simpler depth-averaged representations (Bursik *et al.* 2005). The motion of the material is considered to be gravitationally driven and resisted by both internal and bed friction. The stress boundary conditions are: no stress at the upper free surface and a Coulomb-like friction law imposed at the interface between the material and the basal surface.

The primary factor driving the flow is the component of gravity tangential to the surface, which depends on a local slope computed from the elevation data, hence the criticality of the DEM to the flow computations. The resulting hyperbolic system of equations was solved using a finite-volume scheme with a second-order Godunov solver. Although many real geophysical flows—such as debris flows—are fluidized, this study deals only with granular material that has not been fluidized, such as dome-collapse block and ash flows or rock avalanches initiated by slope instability. The program runs in parallel, using the message passing interface to allow communication between multiple processors, increasing computational power, decreasing computational time and allowing use of large datasets. The algorithm uses local adaptive mesh refinement for shock capturing, and dynamic load balancing for the efficient use of computational resources. Topographic data are included in the simulation through a preprocessing routine in which the digital elevation data are imported. TITAN2D performs flow simulations on a DEM of a desired region, the simulation accuracy being highly dependent on the level of the DEM resolution and quality.

Inputs to the code are the size and location of the initial volume, the internal and bed friction and the DEM. Dalbey *et al.* (2008) presented several methods for characterizing the effect of input data uncertainty on model output. At that time, efficient methods for F representing the uncertainty associated with spatial parameters like terrain elevation were not well understood.

### (f) Bayes linear method

The straightforward way to account for uncertain inputs and stochastic forcing is a Monte Carlo approach—run many simulations and ‘average’ the results in some fashion. If simulations are expensive to run, this approach is not feasible. To circumvent this difficulty, the statistics community has developed the idea of an emulator. In essence, the emulator is a regression surface based on a representative sample of simulations at selected inputs, accompanied by statistical error bounds. Equipped with this surface, output values at new (untested) input values need not be run. Instead, output results can be determined by evaluating the emulator. There are indeed many methods—kriging, metamodels, support vector machines—by which such surrogates may be constructed and there exists a body of literature on the topic (Simpson *et al.* 2001; Clarke *et al.* 2005). One often used emulator is the GAuSsian Process (GASP) emulator, which assumes the regression has the form of a trend plus a Gaussian (Kennedy & O'Hagan 2001; O'Hagan 2006; Bayarri *et al.* 2009; Conti & OHagan 2010). Rougier (2008) in his construction of a multi-variate emulator called the outer product emulator mapped the field output directly by including parametric regression terms on the output index. To construct a GASP emulator, the covariance structure of the Gaussian must be assumed and parameters determined by Bayesian or partially Bayesian methodology. A fully Bayesian determination of the emulator can be costly, especially if the input data are high dimensional. Here, the Bayes linear method (BLM) (Goldstein 1995) to construct an emulator was used. Given prior beliefs (*B*) of mean and variance, the BLM updates these beliefs conditioned on the data (*D*). Note that ‘data’ generally here refers to the output of computationally expensive physics-based simulators. Because only the first two moments of a distribution are determined, the BLM is exact only for Gaussian distributions. As an emulator construction, the BLM update is simpler than a full GASP construction, but the resulting emulator is comparable. Given the prior expectation *E*[*B*] and variance *var*(*B*), the BLM updates are
2.5These update formulae can be derived by minimizing the mean square error (*B*−*a*^{T}*D*)^{2} between *B* and some linear combination of the data. Thus, the BLM update can be viewed as the projection of the set of prior beliefs onto the span of the data.

## 3. Implementation

### (a) Case study I: Galeras Volcano

Galeras Volcano (elevation 4276 m), located in southwestern Colombia (1 ^{°}13.31^{′} N and 77 ^{°}21.68^{′} W), is one of the most active volcanoes in the world (Hurtado & Cortes 1997). Nearly 400 000 people currently live near the volcano; 10 000 of them reside within the zone of high volcanic hazard. Pyroclastic flows pose a major hazard for this population. The current period of activity that began in 2004 (Global Volcanism Program 2012*a*) presents a serious problem for all stakeholders: decision-makers, scientists, public safety officials and the general population. Computational modelling has the potential to provide useful information for hazard assessment and risk mitigation. However, there is a need to evaluate the validity of the modelling and the quality of the DEMs available for use in such modelling.

Galeras is an important volcano for computational flow modelling from both risk management and scientific perspectives (Calvache *et al.* 1997). Forecasts of volcanic explosions using various geophysical tools (Narvaez *et al.* 1997) have occasionally brought public warnings to a high level of alert during the past 20 years. When the alert reaches the highest level, the public are urged to evacuate some local areas; this occurred as recently as in January 2010 (Stefanescu *et al.* 2010*a*; Global Volcanism Program 2012*a*). The worst event at Galeras occurred in 1993, when an eruption killed nine scientists and journalists (Baxter & Gresham 1997).

The topography of the volcano presents a problem for the creation of a good DEM. The irregular morphology on a small scale, with steep slopes, narrow channels, deep gorges and abrupt cliffs, poses problems for the creation of accurate topographic models (Ordoñez Villota & Jentzsch 2000). In addition, the current flow hazard map at Galeras is mainly based on the sparse geological record (Calvache 1990). Dense vegetation, deep erosion, successive deposits of lava and pyroclastic flows hinder the tracing of specific deposits in the field. The diverse effects of this landscape, as reflected in DEMs created by different processes and of different scales, must be examined and quantified to determine the level of confidence that can be placed in model results. Galeras provides a wide range of topographic features that challenge the use of computational flow models.

### (b) Case study II : Mammoth Mountain

Mammoth Mountain is a large, geologically young, composite dome volcano located on the southwestern rim of Long Valley Caldera, CA (Bailey 1989). There are many active hazard issues for Mammoth Mountain, including snow avalanches, rock avalanches and debris flows. In addition, it is intersected by the Mono-Inyo Craters volcanic chain, which is the most active volcanic region in the southwestern USA. If Mono-Inyo-type activity occurs on Mammoth Mountain, then domes may form. These new domes would be growing atop a steep edifice, and therefore could become gravitationally unstable (Hildreth 2004; Global Volcanism Program 2012*b*). Given that block and ash flows occurred at Mammoth Mountain during its older dome growth stage, there is reason to believe that renewed dome formation would result in the activity of block and ash flow. If this is so, then parts of Mammoth Lakes, CA, are at risk from block and ash flows. Mammoth Mountain was used to test the hypothesis that different DEMs result in different model outputs of block and ash flow inundation.

### (c) Model set-up

In quantifying the DEM uncertainties using TITAN2D, a set of parameters was drawn on which to set the bounds of the input domain: internal friction angle, basal friction angle, flow volume, location and DEM. The numerical values for these parameters were chosen to bracket the range of flow volumes and initial locations, and to be representative of the friction angles that have been used by other researchers in their computational models. For the sites used in the study, the surface properties and the rheology are comparable, which is the main reason why the same reasonable parameter values were used for both volcanoes. The internal friction angle has little effect on the output of the flow models (Sheridan *et al.* 2005; Dalbey *et al.* 2008). Many TITAN2D users have chosen values of internal friction that range between 15^{°} and 37^{°} with values between 30^{°} and 35^{°} being the most frequently used (Patra *et al.* 2005; Murcia *et al.* 2010). For the study, an internal friction angle uniformly distributed between 20^{°} and 25^{°} was used.

The value of the basal friction angle has a large effect on flow dynamics in the TITAN2D simulations (Patra *et al.* 2005; Stinton *et al.* 2006). Factors that could affect the choice of basal friction angle include the volume of the flow, the type of pyroclastic flow, the nature of the substrate and the amount of channelization. Murcia *et al.* (2010) listed the basal friction values chosen by TITAN2D users; they range between 5^{°} and 28^{°}, the mean value being about 15^{°}. A basal friction angle uniformly distributed between 15^{°} and 20^{°} was used.

Volumes of pyroclastic flows at stratovolcanoes typically cover a few orders of magnitude. The volume values in this study bracket the range of possible pyroclastic flows for both Mammoth and Galeras. According to Calvache (1990), Galeras Volcano produced five large pyroclastic flow eruptive episodes; a historic eruption in 1866, and prehistoric events in 1100, 2300, 2900 and 4500 yBP. The total deposit volumes of these episodes are in the range . Block and ash flows on Mammoth Mountain might contain –10^{7}) m^{3} of material (Patra *et al.* 2005; Burkett 2007). Thus, the choice of volumes ranges from 1.9×10^{5} to 5×10^{6} m^{3}. The shape of the initial failure region is approximated as a paraboloid of radii , and height . The volume is calculated as For a good match of the volume range, the radius values were uniformly distributed between 25 and 500 m, while the initial height followed the same distribution with values between 10 and 150 m.

Initiation locations were taken from previous mapping of vent sites (Bailey 1989), coupled with knowledge of known weak areas within the volcano as indicated by hydrothermal alteration. Around the centres of the separate initiation locations, different starting positions were uniformly distributed in a circle of radius 200 m. A rectangular area of approximately 40 km^{2} was defined around the vent within the available DEMs as the potential run-out area.

## 4. Results

One of the goals of the analysis was to understand the effect of the spatial structure of available DEMs on hazard maps. Figure 2*c*,*d* shows the correlograms for the ASTER DEM and the TOPSAR DEM, which are the DEMs considered to best represent the real topography for Galeras Volcano and Mammoth Mountain, respectively. It is apparent that data processing resulted in a smoothing and filtering of the TOPSAR DEM, which causes the correlation coefficient to vary smoothly as a function of distance and any two elevation values. Using a distance between two points of 2000 m for the ASTER DEM, the correlation coefficient is 0.6, whereas for the TOPSAR DEM the correlation coefficient is 0.4. This means that elevation values within the ASTER DEM are more highly correlated.

Starting from these premises, the hazard map output for the cases when the DEM is considered to be an input parameter for the TITAN2D model can be explained. Figures 3*a*,*c* and 4*a*,*c* display maps of Galeras and Mammoth, respectively, of the probability that the flow depth will exceed 0.5 m in the next 10 years using methods 1 and 2 to create the terrain realizations.

Figures 3*b*,*d* and 4*b*,*d* show maps at Galeras and Mammoth of the spatially varying lack of confidence in the probability hazard map. The lack of confidence is defined as the computed standard deviation of hazard probability *σ*_{P} divided by the hazard probability, *P*. When calculated by standard means, as was done here, the ratio *σ*_{P}/*P* measures the lack of confidence in the statistic, *P*, owing to insufficient re-sampling of the input parameter space.

Comparison of the figures leads to the conclusion that the difference between hazard map outputs is more pronounced for Galeras Volcano than for Mammoth Mountain. From the lack of confidence figure, it is observed that in both cases the error is concentrated at the flow margins.

For Mammoth Mountain, the differences are less pronounced, but with important differences again concentrated at the edge of the flow. An illustration of how the probabilities vary for method 1 compared with method 2 is shown in figure 5. It was observed in comparing every point where there is a probability of having a flow depth greater than 0.5 m that, the results for Galeras show a greater dispersion than do those for Mammoth Mountain. When the flow is deep, the probability is high and tends to cluster near unity for both mountains. As the probability decreases, dispersion becomes greater for Galeras Volcano. It can be concluded that, as the error map becomes more highly correlated, one should use a more complex method for the creation of topographic realizations such as stochastic method 2. It appears that the spatial autocorrelation of the elevation points influences the hazard map output and a random perturbation of the elevation such as that used in method 1 will not capture this effect.

In previous work (Stefanescu *et al.*, in press), it was concluded that, for moderate and small-sized flows, different representations of the terrain more profoundly affect computation of an accurate flow footprint. For the present contribution, a new set of hazard maps for the case wherein the volume is low was built, with a range between 10^{4} and 5×10^{4} m^{3} and for the high-volume case, using values of 9×10^{6} to 5×10^{7} m^{3}. Since only 517 spatial cells for Mammoth and 872 cells for Galeras were included in the flow footprint for the low-volume case, for any particular cell, the probability that the flow would include that cell tends towards unity (in the case of cells within the starting region) or zero (in the case of any cell outside of the starting region but still within the footprint). Thus, the probability plot is nearly binary, which means that there is a hazard (flow greater than 0.5 m) with either probability approximately 1 or 0 (figure 6*a*,*b*). It can be observed that there is a significant mismatch of prediction between the two methods (left upper corner and right lower corner in both figures) for both volcanoes that can be critical in the case of a hazard or risk assessment. For high-volume flows at Galeras, it was observed that the area of probability greater than zero is much smaller when method 2 (topographic error is spatially correlated) was used, compared with when method 1 was used (no correlation of error between topographic points). This counterintuitive effect needs further study.

The main goal of the study was to explore the effect of DEM uncertainty in constructing a probabilistic hazards map from a geophysical flow model. Quantitative and qualitative analyse were performed for the case wherein the ‘original’ or ‘best’ deterministic DEM (Stefanescu *et al.* 2010*b*) is used as the input parameter for the hierarchical emulator, the output of which is then contrasted with the case, wherein the input is a set of terrain realizations. Hazard maps produced when DEM uncertainty is not included were compared with maps produced when DEM uncertainty is included. Figures 7 and 8 show that, for Galeras, the probability that the flow was deeper than 0.5 m varies considerably from the case of no DEM uncertainty. Hence, the DEM is an important input parameter of which the errors need to be carefully considered in flow modelling, and the effect of the DEM is not diminished by other uncertain parameters or the methodology used. It can be observed that the uncertainty of having flow greater than 0.5 m increases towards the flow edge. For Mammoth Mountain, the DEM uncertainty results in more uncertainty in the flow outline when method 2 is used. One of the causes might be that the flow propagates a shorter distance than that in the original DEM. At Galeras, where the autocorrelation is high, the uncertainty in the flow outline increases when correlated DEM error is taken into account. In this case, the uncertainty in flow outline increases, suggesting that perturbing the DEM is more important as autocorrelation increases.

## 5. Conclusions

Computer models of hazardous phenomena, such as floods, hurricanes and avalanches, are expensive to run, and each run produces an enormous amount of data. For example, a flood model output may consist of water depth and velocity at every point in a large grid, at every time step. Furthermore, these models often require specification of several parameters that may not be well characterized, and initial and boundary data that are likewise poorly specified. For the first time, in this contribution, a process for computing a hazard map owing to a geophysical flow with both uncertain model input parameters and an uncertain elevation map has been described. Uncertainty in the elevation map has been addressed using two methods for creating an ensemble of map realizations with the same error structure as an original elevation model. In one method, the errors in the model are assumed to be spatially uncorrelated, whereas, in the other method, the autocorrelation structure of the error was used to produce the realizations. Once the elevation realizations are produced, the computational model is run over the ensemble of elevation surfaces and outputs are appropriately combined. The results suggest that it is critical to consider the error autocorrelation structure in the DEM to properly incorporate DEM error in the entire error model and the resulting probabilistic hazards map.

For the two test sites used in this study, one of the main differences is the texture of the terrain surface: the digital representation of the surface of Galeras Volcano is quantifiably rougher than is the representation of Mammoth Mountain. An important conclusion for researchers is that, based on the surface roughness of the area of study, different methods to assess the uncertainty caused by the DEM in the flow model can be implemented. In the case of a smooth terrain, for a fast and less expensive computational implementation, a simple ‘random’ perturbation method (method 1) yields results similar to those using the stochastic method (method 2). For rough surfaces, the method for creation of separate terrain realizations should include some characterization of the autocorrelation structure in the DEM error. A method based on the error map and an unconditional stochastic simulation as presented in this paper is a good option.

## Acknowledgements

This work was supported by NASA grant NNX08AF75G. The work and opinions expressed herein are those of the authors alone and do not reflect the opinion of NASA. We are grateful to JPL for the construction and distribution of the TOPSAR dataset.

- Received December 6, 2011.
- Accepted January 12, 2012.

- This journal is © 2012 The Royal Society