The sensitivity of a volcanic flow model to digital elevation model accuracy: experiments with digitised map contours and interferometric SAR at Ruapehu and Taranaki volcanoes, New Zealand

Journal of Volcanology and Geothermal Research 119 (2002) 89^105 www.elsevier.com/locate/jvolgeores

The sensitivity of a volcanic £ow model to digital elevation model accuracy: experiments with digitised map contours and interferometric SAR at Ruapehu and Taranaki volcanoes, New Zealand N.F. Stevens a;
, V. Manville b , D.W. Heron a a

Institute of Geological and Nuclear Sciences, Grace¢eld Research Centre, P.O. Box 30368, Lower Hutt, New Zealand b Institute of Geological and Nuclear Sciences, Wairakei Research Centre, Private Bag 2000, Taupo, New Zealand Received 4 September 2001; accepted 25 March 2002

Abstract A growing trend in the field of volcanic hazard assessment is the use of computer models of a variety of flows to predict potential areas of devastation. These can be compared against historic and geologic evidence of past events for model calibration, or used to construct hazard zone maps for mitigation and planning purposes. The accuracy of these computer models depends on two factors, the nature and veracity of the flow model itself, and the accuracy of the topographic data set over which it is run. All digital elevation models (DEMs) contain innate errors. The nature of these depends on the accuracy of the original measurements of the terrain, and on the method used to build the DEM. In this paper we investigate the effect that these errors have on the performance of a volcanic flow model designed to delineate areas at risk from lahar inundation. The model was run over two DEMs of southern Ruapehu volcano derived from (1) digitised 1:50 000 topographic maps, and (2) airborne C-band synthetic aperture radar (SAR) interferometry obtained using the NASA AIRSAR system. On steep slopes of V4‡ or more, drainage channels are more likely to be incised deeply, and flow paths predicted by the model are generally in agreement for both DEMs despite the differing nature of the source data. Over shallow slopes (V 6 4‡), where channels are less deep and are more likely to meander, problems were encountered with flow path prediction in both DEMs due to interpolation errors between contours, and due to forestry. The predicted lateral and longitudinal extent of deposit inundation was also sensitive to the type of DEM used, most likely in response to the differing degrees of surface texture preserved in the DEMs. A technique to refine contour-derived DEMs and reduce the error in predicted flow paths was tested to improve the reliability of the modelled flow path predictions. In areas where high-resolution topographic maps are unavailable, forthcoming topographic measurements acquired by a single-pass space-borne instrument, the NASA Shuttle Radar Topography Mission (SRTM) are likely to prove invaluable. 8 2002 Elsevier Science B.V. All rights reserved. Keywords: digital elevation model; LAHARZ £ow model; debris avalanche; lahar; topographic map; interferometric synthetic aperture radar; Ruapehu; Taranaki/Mt Egmont

* Corresponding author. E-mail address: [email protected] (N.F. Stevens).

0377-0273 / 02 / $ ^ see front matter 8 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 7 3 ( 0 2 ) 0 0 3 0 7 - 4

VOLGEO 2507 9-11-02

90

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

1. Introduction Gravity-driven £ows constitute a signi¢cant hazard at numerous volcanoes world-wide, including such phenomena as lava £ows, blockand-ash £ows, and lahars and debris avalanches. Traditional methods of assessing the hazard zones associated with these events are based on reviews of historical records and ¢eldwork to identify the limits of their deposits in the geological record. Predictions of future hazard zones are then based on interpolation and extrapolation of known data, perhaps supplemented by calibrated £owrouting models (Costa, 1997; Laenen and Hansen, 1988; Macedonio and Pareschi, 1992), and the researchers’ own experience. However, these techniques are less useful at volcanoes where little data on past activity are available, so that predictive models which calculate hazard zones from either the physics of the £ow (Malin and Sheridan, 1982), or empirical relationships between the size of the £ow and the impacted area (Iverson et al., 1998), become more important. The gravitational collapse of the upper parts of volcanic edi¢ces to form massive debris avalanches and long run-out lahars (Glicken, 1982, 1998; Crandell et al., 1984; Scott and Vallance, 1995) is one of the most serious of volcanic hazards due to their extreme mobility (Carrasco-Nunez et al., 1993; Vallance and Scott, 1997). A method of delineating areas at volcanoes vulnerable to debris avalanche hazards is computer £ow modelling (e.g. Iverson et al., 1998; Schilling, 1998). This is achieved by calculating the predicted £ow path and deposit inundation over a representation of the local topography in the form of a digital elevation model (DEM). Reasonably accurate, up-to-date topographic data are therefore invaluable for this type of risk assessment. However, all DEMs contain innate errors (e.g. Shearer, 1990; Wechsler, 2000), and to date, little work has been done to assess the sensitivity of the predictions of volcanic £ow models against DEM quality. In this paper, we investigate the sensitivity of a volcanic £ow model, LAHARZ, developed at the United States Geological Survey to model debris avalanches and long run-out lahars (Iverson et al.,

1998; Schilling, 1998), to DEM quality. We use DEMs from two diverse sources as input: (1) a DEM generated from the national 1:50 000 topographic map series, and (2) data acquired during an aerial interferometric radar campaign (NASA TOPSAR). We identify disparities in the £ow model output, as a result of running the model over the di¡erent DEMs, and discover the causes by identifying particular types of error in the DEMs. We develop methods to double-check and improve the modelled £ow paths. A rationale for this study is to ascertain whether improved topographic mapping is necessary for realistic computerised risk assessment at other volcanoes for which 1:50 000 maps are the only topographic data available (e.g. Taranaki/Mt Egmont, New Zealand). The suitability of a forthcoming topographic data set of almost world-wide coverage, namely the NASA Space Shuttle Radar Topography Mission (SRTM) data (van Zyl, 2001; Werner, 2001), is also discussed and evaluated theoretically using resampled TOPSAR data.

2. The study area New Zealand contains a number of active andesitic stratovolcanoes, including Ruapehu, Ngauruhoe and Tongariro in the Tongariro Volcanic Centre at the southern end of the Taupo Volcanic Zone (Fig. 1), and Taranaki/Mt Egmont 130 km to the west. The composite cone of Ruapehu is the highest peak in the North Island of New Zealand, comprising a 110-km3 volcanic edi¢ce (Hackett and Houghton, 1989), surrounded by a radial ring plain of similar volume (Cronin and Neall, 1997; Cronin et al., 1996; Lecointre et al., 1998; Palmer, 1991; Palmer et al., 1993) that includes at least one debris avalanche deposit (Palmer and Neall, 1989). The summit area supports a number of perennial snow¢elds and glaciers while the active crater hosts a hot acidic lake (Christenson and Wood, 1993), making the generation of lahars a particular feature of volcanic activity at Ruapehu. Prehistoric and historic lahars at Ruapehu have been triggered by a variety of mechanisms including explosive ejection of Crater Lake water (Cronin et al., 1997a; Healy

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

91

Fig. 1. Location map showing clockwise from top left, the North Island of New Zealand, the Tongariro Volcanic Centre and adjacent infrastructure (Ruapehu shown by light shading, Tongariro volcano by darker shades), and a shaded relief map of Ruapehu, Tongariro, and the cone of Ngauruhoe overlain by a drainage network derived from aerial photography. Ruapehu volcano is bordered to the west and east by major normal fault scarps and mountain ranges, and to the north by Tongariro. These surrounding topographic features in£uence the drainage pattern on the ring plain. Rivers on the southeastern £anks curve westwards with decreasing elevation rather than following radial paths. The dashed lines mark the edges of the swath of 1996 NASA TOPSAR data used in this study. The white rectangle shows the extent of the study area.

VOLGEO 2507 9-11-02

92

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

et al., 1978; Nairn et al., 1979), remobilisation of tephra deposits from the cone (Manville et al., 2000), and the failure of tephra barriers at the Crater Lake outlet (O’Shea, 1954). Ruapehu, and the rest of the Tongariro Volcanic Centre, is unusual by New Zealand standards because it is covered by two sets of digital topographic data obtained by di¡erent mapping and data processing methods, and are of di¡erent levels of accuracy. Firstly, the whole of New Zealand is mapped at 1:50 000 for the national topographic map series (e.g. Land Information New Zealand, 1987, 1996a,b, 1997, 1998). These data were released recently in digitised format, and are used to derive DEMs by various interpolation methods. Secondly, the area was mapped topographically in November 1996 by the NASA TOPSAR interferometric synthetic aperture radar (SAR) system, which was carried onboard a DC-8 aircraft. These data are referred to here as the TOPSAR DEM. The latter has a superior accuracy, but is limited in spatial extent, being obtained in 10-km swaths. The southern Ruapehu volcano area was chosen as a ¢eld site because both Land Information New Zealand (LINZ) and TOPSAR DEMs cover a signi¢cant length of the major river catchments, radar shadowing e¡ects are minimal, and because the area is not the focus of current lahar activity. The volcano consists of two distinct, contrasting types of terrain, (1) the £anks of the cone, which are characterised by relatively steep, deeply incised channels; and (2) the shallower slopes (V 6 4‡) of the surrounding volcaniclastic ring plain (Palmer et al., 1993), where drainage channels are less distinct and more sinuous, and where most of the local infrastructure is sited.

3. The LAHARZ model The LAHARZ model was developed to delineate potential inundation zones, using a DEM, user-speci¢ed debris volumes, and the ARC/INFO0 Geographic Information System (Schilling, 1998). LAHARZ was designed to provide a fast, objective, automated and reproducible method of risk assessment (Iverson et al., 1998).

The processing sequence is described brie£y here to aid later discussion. The DEM is ¢rst processed to ¢ll sinks, or depressions. The surface hydrology functions in the ARC/INFO0 GRID software module are then used to derive the £ow paths and £ow direction across the DEM. These are identi¢ed by a user-de¢ned threshold number ^ if the number of cells ‘£owing’ into an individual cell exceeds the threshold number, then the cell is de¢ned as being part of one of the £ow paths. The hazard zone is calculated using empirical statistical relationships derived by Iverson et al. (1998) between lahar volume, and crosssectional and planimetric area of lahar inundation. Inundation is calculated at each cell along a £ow path, from the edge of a calculated proximal hazard zone (Malin and Sheridan, 1982) by ¢lling the corresponding valley cross-section for that cell to the level required to satisfy the empirical relationship. A running tally is kept of the volume of material ‘stored’ in the ¢lled cells so that the lahar only propagates until it runs out of material. The DEM and four user-speci¢ed volumes are used to delineate degrees of lahar inundation. We have speci¢ed volumes in the simulations in this study ranging from 1.0U107 to 2.5U108 m3 . These were chosen to re£ect a reasonable range of potential activity, and it is stressed here that the impetus for our study was not general risk assessment for the area, which is the subject of on-going research. The main focus was on evaluation of model sensitivity to DEM quality, and although lahars at Ruapehu with volumes outside of this range do occur (Cronin et al., 1997b; Palmer and Neall, 1989), these values are adequate for our evaluation purposes. It is beyond the scope of this paper to discuss the general performance of LAHARZ relative to real £ow events and physical processes.

4. The sources of digital elevation data, LINZ and TOPSAR Of the two sets of digital elevation data available for use in this study, the LINZ data represent the publicly available, low-cost option in New

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

93

Fig. 2. (a) Shaded relief image of the LINZ DEM. (b) Close-up of the surface of the LINZ DEM at the distal end of the Tahatekapua lava £ow. The £at terraces in the shallow area to the south of the lava £ow are caused by interpolation between points at the same elevation. (c) Close-up of the TOPSAR DEM in the same area. Elevation measurements at 10-m intervals preserve surface detail, including sinuous channels, ¢eld boundaries (indicated by black arrow) and the railway cutting (marked by white arrow). Although the DEMs are processed to a similar spatial resolution, the ¢ne detail is best preserved in the TOPSAR data.

Zealand at the time of writing. The LINZ data set consists of digitised 20-m contours and spot heights, derived from 1:50 000 topographic maps (Land Information New Zealand, 1987, 1996a,b, 1997, 1998), and also includes vector layers of the rivers and other features included in the topographic maps (Fig. 1). A DEM of 20-m horizontal resolution was generated from the contour and spot height data, and is referred to here as the LINZ DEM (Fig. 2a,b). The TOPSAR DEM was commissioned during the November 1996 NASA PACRIM mission, and was obtained using airborne interferometric SAR (Zebker et al., 1992). Topography is derived from phase shifts in simultaneous SAR acquisitions, measured by two spatially separated anten-

nae mounted on the aircraft fuselage (Fig. 3a). The technique is already described in detail elsewhere (e.g. Rowland et al., 1999, Appendix I), and we will not elaborate further here. The resulting DEM was processed to a horizontal resolution of 10 m (Fig. 2c) and georeferenced. The di¡erent methods used to generate these DEMs means their overall accuracy will be di¡erent, and they will contain errors arising from diverse sources. The ultimate accuracy of DEMs derived from digitised contours, as in the LINZ DEM, depends on (a) the character of the terrain, (b) the contour interval, (c) the scale of the topographic map, and (d) the method used to interpolate between contours (Robinson, 1994). Previous investigations into the accuracy of the LINZ

VOLGEO 2507 9-11-02

94

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

Fig. 3. Single-pass interferometric SAR data are acquired by transmitting radar energy towards the Earth’s surface from an onboard instrument (solid arrow), and by measuring the re£ected signal (dashed arrows) at two onboard spatially separated antenna, labelled ‘1’ and ‘2’ here. (a) The NASA AIRSAR system has antennae mounted on the fuselage of a DC-8 aircraft, separated by 2.5 m. (b) To obtain suitable antennae separation from the orbiting altitude of the Shuttle, the second antenna was placed at the end of a 60-m boom, which extended from the payload bay.

DEM (Barringer and Lilburne, 1997), as well as assessments of DEMs derived from digitised map contours in general (e.g. Shearer, 1990), estimate that this DEM has a root mean squared (rms) vertical accuracy of V8 m, although this would vary with slope and local terrain. However, LINZ estimates the horizontal accuracy of the 1:50 000 topographic maps to be no better than N 22 m and the vertical accuracy to be no better than N 10 m (Land Information New Zealand, 2000), so the actual vertical DEM accuracy will probably fall short of this V8-m rms estimate. It is emphasised here that decreasing the horizontal processing resolution is unlikely to improve the accuracy of a DEM derived from map contours, if this has been chosen appropriately to begin with (Shearer, 1990). ARC/INFO0 provides two techniques for the creation of DEMs from digitised contours. In earlier versions of the software, DEMs were created from a triangular irregular network (TIN) generated from contour and spot height data. DEMs derived from TINs often contain artefacts of zero slope, particularly in areas of sinuous widely

spaced contours. These occur where the three points that de¢ne a TIN triangle are derived from the same contour (Fig. 4). This results in £at terraces along the position of the original contour lines that are referred to as ‘£at spots’ by Robinson (1994). ARC/INFO0 version 7 provided the Topogrid interpolation method based on the ANUDEM program developed by Hutchinson (1988, 1989) and designed for the creation of hydrologically correct DEMs. Topogrid uses iterative ¢nite di¡erence interpolation and is essentially a discretised thin plate spline technique (Wahba, 1990). Topogrid, however, also produces £attening of the elevation model as it crosses the contour (Fig. 2b). In contrast, the TOPSAR DEM is derived from a direct measurement of elevation within every 10m pixel where a coherent radar return is obtained, and has an estimated rms vertical accuracy of 1^3 m in mountainous terrain (Madsen et al., 1995). The TOPSAR DEM is thus not a¡ected by the interpolation errors a¡ecting the LINZ DEM (Fig. 2c), and high-resolution surface texture is recorded. The data are however a¡ected in local-

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

95

5. Comparison between model output with the LINZ and TOPSAR DEMs as input data The LAHARZ model was run using the LINZ and TOPSAR DEMs of southern Ruapehu as inputs. In each case this was a two-fold process, examining ¢rst the £ow path prediction, and second, the calculated lahar inundation. 5.1. The modelled £ow paths

Fig. 4. (a) Schematic diagram showing how £at spots occur in DEMs derived from topographic maps. Three contours (thick curved lines) are digitised and are stored as vectors using a series of digitisation points or nodes (grey dots). When a TIN is generated between the nodes to represent the surface of the landscape, a number of £at spots are generated where nodes of equal height are joined (grey triangles), forming horizontal terraces along the position of the contours. These may be interpreted as hydrological £ow paths during automated £ow path calculation.

ised areas by radar shadowing. This occurs when the radar return to the instrument is weak, generally in slopes facing away from the sensor or in areas obscured from the line-of-sight of the radar (e.g. at the bottom of deep, narrow channels). When small areas are a¡ected, it is possible to interpolate across these with reasonable accuracy, providing the terrain is fairly uniform. In this study, Delaunay triangulation was used to interpolate across areas a¡ected by radar shadow of less than ¢ve pixels across. A problem associated with both DEMs occurs in forested regions, where the height of the canopy is measured, rather than that of the forest £oor. This is apparent in the shaded relief image of the TOPSAR DEM in Fig. 2c, where the linear ¢eld boundaries can be seen clearly in the southern- and eastern-most parts of the image. Compensations can be made in contours derived from stereo photographs in topographic maps, but not always (e.g. Stevens et al., 1999), and undulations in the canopy height cannot be anticipated.

Flow paths were calculated across the LINZ and TOPSAR DEMs and are superimposed for comparison in Fig. 5. On the steep slopes of the volcano £anks (the top half of Fig. 5), the modelled £ow paths were in close agreement. This is not altogether surprising, because the channels in this area tend to be deeply incised and wide, and the contours of the LINZ data are closely spaced, meaning the interpolation between contours is most likely to be successful. In one instance (marked with an ‘s’ in Fig. 5), a section of the £ow path was not identi¢ed in the TOPSAR data. This was found to be due to a remaining area of radar shadowing, which had not been removed by interpolation, and whose masking value formed a ‘hole’ in the DEM. Over the lower gradient ring plain (the lower half of Fig. 5), there was signi¢cant disagreement between the modelled £ow paths. We investigated the cause of these by close examination of the DEMs, and found they are caused directly by the aforementioned £aws in the DEMs. Where the modelled £ow paths run parallel to and within a few 100s of metres of the actual channel position overall, and rejoin the same £ow path, this does not compromise the direction of the predicted lahar path signi¢cantly. The nature of the LINZ DEM (Fig. 2b) makes it improbable that every undulation of a drainage channel would be recreated, especially in the ring plain, where contours can be separated spatially by 750 m and more, and where drainage channels are increasingly shallow and sinuous. Planar facets created during the interpolation process can create unusual modelled drainage patterns (Fig. 5, marked ‘t’). Of greater concern were instances where £ow

VOLGEO 2507 9-11-02

96

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

Fig. 5. Flow paths calculated with LAHARZ using the LINZ DEM (black lines) and the TOPSAR DEM (white lines) for southern Ruapehu. The position of the Tahatekapua lava £ow, two large fault scarps and the railway cutting are shown in light grey for orientation. Di¡erences in predicted £ow paths are due to, C ^ erroneous perpendicular £ow path created by interpolation of LINZ data; V ^ erroneous perpendicular £ow path created by forest canopy in TOPSAR data; S ^ break in £ow path due to radar shadowing in TOPSAR data; and T ^ parallel drainage resulting from planar facets in the DEM produced during interpolation.

paths changed direction abruptly by 90‡ and ran perpendicular to and across the drainage network derived from the other DEM (Fig. 5, marked by ‘c’ and ‘v’). To double-check the validity of the

routes of the modelled £ow paths, we compared Fig. 5 with the vector layer of the streams in the area (Land Information New Zealand, 1987, 1996a,b, 1997, 1998). These were mapped origi-

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

97

Fig. 6. Predicted lahar inundation using the same starting position in (a) the LINZ DEM, and (b) the TOPSAR DEM, plotted over shaded relief images of the respective DEMs at the same scale. Because a £at spot was identi¢ed incorrectly in the LINZ DEM as a £ow path, LAHARZ delineates di¡erent areas at risk. Arrows indicate the major fault scarps, and the Tahatekapua lava £ow is in the centre of each DEM. ‘X’ indicates areas in the LINZ DEM where the lahar inundation has spread unreasonably far laterally, coinciding with £at spots in the DEM. ‘f’ denotes an area of forest in the TOPSAR DEM, which is incorrectly interpreted as a high relief obstacle to the lahar deposition.

nally from aerial photographs and can be considered representative of the actual drainage pattern, although this will be discussed further, later. The TOPSAR-derived £ow paths are most similar to the rivers in the vector layer, and also most closely follow the sinuousness of the actual drainage network. All but one of the erroneous perpendicular £ow paths were created using the LINZ DEM (marked ‘c’ in Fig. 5). These all cut across two or more channels, and joined a parallel channel kilometres away from the original stream course. This would mean a modelled lahar would jump incorrectly to another part of the drainage network, and the predicted area of inundation downstream would be signi¢cantly di¡erent (Fig. 6). Comparison of the erroneous perpendicular LINZ-derived £ow paths with a vector layer of the original 1:50 000 con-

tours shows that the errors coincide with and follow a particularly sinuous contour. The horizontal terraces caused by £at spots in the LINZ DEM that arise from incorrect contour interpolation are therefore the most likely cause of this e¡ect. Less straightforward to explain was the one instance of this e¡ect occurring in the TOPSARderived £ow paths (marked ‘v’ in Fig. 5). Comparison with the 1:50 000 LINZ topographic map showed that the erroneous perpendicular £ow path coincides with a ¢eld boundary, which is visible in a shaded relief image of the TOPSAR DEM as the edge of a forest. According to the map, the river £ows into the forest, but in the TOPSAR data, the forest canopy obscures the river channel. The e¡ect was not found at the numerous other forestry boundaries present in the area.

VOLGEO 2507 9-11-02

98

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

Fig. 7. Predicted lahar inundation from the same starting point in (a) the LINZ DEM, and (b) the TOPSAR DEM. In this case, the modelled £ow paths are in agreement, but LAHARZ delineates slightly di¡erent areas at risk because of the way the deposits are apportioned upstream. The arrows and ‘X’ have the same assignation as in Fig. 6.

5.2. The modelled lahar inundation

6. Evaluating £ow paths at Taranaki from LINZ digital elevation data

Lahar inundation was modelled to the west of the Tahatekapua lava £ow over both DEMs using £ow paths that had been located correctly in both DEMs (Fig. 7). The LINZ-derived inundation is spread wider laterally (Figs. 6a and 7a, marked ‘x’), so for all four user-speci¢ed lahar volumes, the lahar does not reach the same length as the TOPSAR-derived inundation. This is an e¡ect of the horizontal terraces along LINZ contours, which allow the modelled deposits to spread further laterally and reduce the lahar volume available for downstream propagation. The modelled deposition is also ‘blocked’ by the forest canopy in the TOPSAR DEM (Fig. 6b, marked ‘f’). Both DEMs provide a similar delineation of the lahar deposit, although in general, the TOPSAR data are better for ¢ne topographic detail (e.g. observed interaction between the lava margin, a fault scarp and the modelled lahar margin in Fig. 6b).

Part of the impetus for this study was to ascertain whether volcanoes, where topography is mapped to 1:50 000 at best, need to be mapped in ¢ner detail to achieve reliable modelling for risk assessment. The Taranaki volcanic centre (Fig. 8), 130 km to the west of Ruapehu, consists of a group of four volcanoes including the symmetrical young cone of Taranaki/Mt Egmont (Neall et al., 1986). The current edi¢ce is less than 7000 yr old and is surrounded by a radial ring plain built up of at least 10 debris avalanches and other interbedded volcaniclastic deposits during the volcano’s 130-ka history (Palmer and Neall, 1991; Ui et al., 1986). Although the area is covered by the LINZ 1:50 000 topographic map series, TOPSAR data, or data equivalent in resolution and accuracy, were not acquired at the time of writing, and in this re-

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

99

Fig. 8. Shaded relief model of Taranaki/Mt Egmont overlain by a drainage network derived from aerial photography. Note the radial, almost symmetrical pattern of drainage at Taranaki relative to the more complex drainage system at Ruapehu in Fig. 1. White triangle in inset shows location of Taranaki.

spect Taranaki is similar to many volcanoes around the world. Flow paths were simulated across a DEM of Taranaki derived from LINZ data using LAHARZ, and were compared with the LINZ river vector layer for validation. Places where these differed enough to cause a simulated lahar to £ow into a di¡erent part of the catchment were noted, and were compared with a slope map derived

from the DEM. Of the 24 instances noted, all but three occurred on slopes of 3‡ or less, and none occurred on slopes s 6‡. The LINZ DEM and vector data were examined to determine the cause of the di¡erences between the modelled £ow paths and the actual drainage channels. In three cases it was not possible to determine the cause. Of the remaining 21 instances, nine were caused by £at spots either being interpreted as channels,

VOLGEO 2507 9-11-02

100

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

or diverting the £ow path enough to a¡ect down slope pixels. Twelve di¡erences appear to have been caused by anomalies in the DEM during interpolation between contours. Flat spots were the main cause of erroneous £ow paths at Ruapehu, but this was obviously not the case at Taranaki. Where erroneous £ow paths were calculated for Taranaki, these connected generally to adjacent streams, and did not cross two or more adjacent, parallel £ow paths, as at Ruapehu (Fig. 5, marked ‘c’). This is most likely because the drainage network at Ruapehu is more complex relative to the terrain than at Taranaki (Figs. 1 and 8). Southern Ruapehu has an unusual drainage pattern, with rivers spiralling westwards across its £anks rather than following the steepest route to the ring plain, so that river channels are not perpendicular to the contour lines, which may cause errors due to £at spots to dominate.

7. Preserving £ow paths in contour-derived DEMs Interpolation problems in DEMs derived from topographic maps can be reduced by a number of methods. Although topographic maps at ¢ner resolutions may aid interpolation between contours, £at spots will still a¡ect sinuous contours (e.g. Stevens et al., 1999). Flat spots can be removed manually at the TIN stage of processing, if used (e.g. Robinson, 1994; Stevens et al., 1999), by swapping the vertices of the triangles (Fig. 9A), but this process is prohibitively timeconsuming in large DEMs. Similarly, spot heights or intermediate contours can be inserted before the triangulation stage, but this is also time-consuming. Streams, ridges, and other linear topographic features can be used to create breaklines during the TIN construction, but of these, only the stream data are readily available (Fig. 9B). The ARC/INFO0 Topogrid module includes modi¢cations to allow the ¢tted DEM to follow abrupt changes in terrain. Topogrid uses the contours to generate a network of streams and ridges. Stream data, if available, may be used to constrain the interpolation process. The

Fig. 9. (A) Diagram showing the e¡ect of swapping the vertices between vector nodes to remove £at spots in DEMs (see Fig. 4). This may be achieved manually, but it is labour-intensive. (B) Vector data delineating the position of the stream (denoted by the grey dashed line) are used as ‘breaklines’ to constrain triangulation over the stream channel. Although £at spots are removed along the stream, they remain in other parts of the DEM that are unconstrained by breaklines.

LINZ DEM of Ruapehu was reconstructed using the LINZ stream vectors as breaklines. This reduced the £attening close to the stream channel, but did not eliminate it in the inter-channel areas (Fig. 9B). LAHARZ was run across the breakline-constrained DEM, and the modelled £ow paths agreed well with both the LINZ stream vector layer, and the TOPSAR-modelled channels. The only di¡erences noted were caused by

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

irregular interpolation anomalies, which are more problematic to deal with. These did not cause predicted £ow to join other catchment areas. There were no signi¢cant di¡erences in the predicted LINZ lahar inundation compared with prior attempts when stream vector data were not used as breaklines.

8. Future prospects using SRTM data During February 2000, the NASA Shuttle Radar Topography Mission was completed successfully (SRTM Project O⁄ce, 2001; van Zyl, 2001; Werner, 2001). The Space Shuttle carried a singlepass interferometric SAR, with antenna separation achieved using a 60-m boom extending from the payload bay (Fig. 3b). Topographic mapping in 50-km swaths of all landmass between latitudes 60‡N and 56‡S was achieved using Cband radar, apart from some areas within the USA, with a less extensive, high-resolution coverage being achieved simultaneously using X-band radar. The data are being processed currently to two levels. Level 1 data sets will be processed to 3 arc second resolution (V90 m) and will be available freely. Level 2 data will be processed to 30-m resolution, with distribution of the C-band data likely at the discretion of the US Department of Defense (NASA/NIMA Memorandum of Understanding, summarised, 2000), although this appears to be under discussion at the time of writing. This data set may provide a future means to circumvent the problems associated with deriving £ow paths over DEMs derived from topographic maps. To evaluate the utility of SRTM topographic data for volcanic £ow risk assessment, we resampled the TOPSAR data to 30- and 90m grid resolutions using bilinear interpolation. The ¢rst observed e¡ect was that increasing the pixel size decreased the ¢le size, reducing processing time dramatically. Assuming a 2-m vertical rms error in the original TOPSAR data (Madsen et al., 1995), the rms error of the resampled DEMs was calculated by subtracting the original TOPSAR DEM from the resampled DEMs and using the relation:

rmserror

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N uX u ðh3h Þ2 T u t 1 ¼ N31

101

ð1Þ

where N is the number of sample points (in this case all valid measurements of elevation in the TOPSAR DEM), h is the pixel elevation of the resampled DEM, and hT is the pixel elevation of the TOPSAR DEM. The rms vertical errors are 3.4 and 8.2 m, in the 30- and 90-m resampled DEMs respectively. The vertical accuracy of the resampled TOPSAR data is similar to the expected relative vertical accuracy of the SRTM data ( 6 6 m and 6 10 m in the X- and C-band level 2 data respectively (Deutsches Zentrum fu«r Luft- und Raumfahrt, 1999)). This evaluation of the suitability of the SRTM DEMs is focussed towards feature preservation in grids of 30- and 90-m pixel size at this expected relative vertical accuracy, and the e¡ect of this on LAHARZ. The £ow paths simulated by LAHARZ across the DEMs were compared. Flow paths coinciding with the steeper slopes of the volcano £anks were similar despite resampling, although their sinuousness was generally reduced with increasing pixel size. One case of signi¢cant disparity between £ow paths was noted at the neck of the Tahatekapua lava £ow, where the gradient is low, and the £ow path in the 90-m DEM is predicted to £ow south instead of west. In this case, it is anticipated that surface detail that constrained the £ow paths in the ¢ne resolution DEMs has been lost. Errors such as these, although rare, would therefore need to be £agged by comparison with stream vector data as described before. On the low gradient ring plain, there were four minor disagreements between predicted £ow paths that linked tributaries incorrectly to the main channel, but none of these routed the lahar path into other parts of the catchment. The forest boundary that was incorrectly £agged as a £ow in the TOPSAR data (Fig. 5, marked ‘v’) was also incorrectly interpreted as such in the resampled DEMs. Three volumetrically identical lahars were simulated across the resampled DEMs from the same starting point (Fig. 10). The simulations with vol-

VOLGEO 2507 9-11-02

102

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

Fig. 10. The e¡ect of resampling the TOPSAR DEM (a) to coarser resolutions of 30 m and 90 m (b and c), to simulate SRTM data has little e¡ect on the lahar inundation predicted by the LAHARZ model, other than that the smaller lahars are slightly shorter in length.

umes 2.5U108 m3 and 108 m3 are reasonably similar in lateral and longitudinal extent across all DEMs, but the 107 -m3 volume lahars fall short in length by 1^2 km in the 30- and 90-m resolution DEMs respectively (Fig. 10). This indicates that the deposits are allocated di¡erently per cell upstream, possibly due to the smoothing e¡ect on the terrain of the resampled pixel resolution.

9. Discussion That the majority of the errors in £ow path calculation during this study occurred on the ring plain is a concern, because at many volcanoes, this coincides with the highest density of local infrastructure and settlement (e.g. Fig. 1). The stream vector layer provided a reasonable method for double-checking the modelled £ow paths, and in the case of the contour-derived DEMs, the stream vectors provided a method of

correcting modelled £ow paths. Georeferenced aerial photographs or remotely sensed images could equally be used. A di⁄culty exists when stream courses are not included in the stream vector layer. A clear example of this is evident in the LINZ data of Taranaki (Fig. 8), where the density of mapped streams becomes markedly less on the upper slopes of Taranaki. This coincides with the border of a National Park and an abrupt change in land cover from pasture to native forest (Pairman et al., 1997). Some of the streams were thus not distinguishable in the original aerial photographs in areas covered by dense vegetation. Both approaches to mapping topography examined in this study are a¡ected adversely by the presence of commercial and native forestry. This is most apparent in the TOPSAR DEM because the ¢ne data resolution highlights the edges of commercial forests in southern Ruapehu (Fig. 2c), and this will also be a problem in the forth-

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

coming SRTM data. Forest canopies are not necessarily accounted for when deriving the contours of topographic maps from aerial photographs either (e.g. Stevens et al., 1999). Even if the tree canopy were accounted for, it would be di⁄cult to correct precisely for the variability in forest canopy height, especially in river channels (cf. Downs and Simon, 2001). Methods of topographic mapping using instruments sensitive enough to measure the ground surface beneath the biomass, e.g. laser altimetry, could solve this problem (e.g. Harding et al., 2001). Although we did not address the issue directly during this study using a multi-temporal set of DEMs, we note here that changes in topography at volcanoes, due to deposition and erosion following large-scale lahars and other volcanic processes, are another source of error in digital elevation data. Existing £ow paths may change signi¢cantly as a result. Repeat surveys of volcanoes at a reasonable frequency to record major topographic alteration are recommended.

10. Conclusions Now that we have an understanding of the effect of DEM errors on the volcanic £ow model, LAHARZ, we are in a stronger position to interpret its output. DEMs derived from interpolated contours, as well as those measured in a spatially continuous array, e.g. by interferometric radar, contain innate errors. We have investigated the e¡ect of these errors on £ow path computation and the delineation of volcanic £ow risk, and ¢nd that LAHARZ is most sensitive to DEM errors in shallowly sloping areas (V 6 4‡). This ¢nding also has important implications for similar models used for non-volcanic applications, e.g. landslides and hydrologic transport of pollutants. The quality of the £ow path computation can be double-checked readily to ensure realistic risk delineation, using a map of the actual stream channels, digitised from topographic maps or derived from remotely sensed data. A vector map of the stream channels can also be used to improve interpolation when generating a DEM from digitised contours.

103

The modelling of £ow paths on volcanoes with a strongly radial drainage pattern, such as Taranaki, is less sensitive to errors in the DEM than in volcanoes with complex drainage paths, such as Ruapehu. To a lesser extent, the quality and resolution of the DEM a¡ects the calculated areas of inundation using the LAHARZ model. Finally, we have shown that digital elevation data at the spatial resolution and relative accuracy of the NASA SRTM data will be a valuable and reasonably reliable asset for future risk assessment. SRTM provides a world-wide baseline data set for the year 2001. A future issue of concern will be ensuring these measurements are updated following signi¢cant topographic alteration due to volcanic activity, and other large-scale processes. An airborne system similar to the NASA interferometric SAR, TOPSAR, would be ideal for this purpose.

Acknowledgements This project was supported by funding from the New Zealand Foundation for Research in Science and Technology (FRST). We are grateful to Ellen O’Leary and team, NASA Jet Propulsion Laboratory, for the timely provision of the TOPSAR data. The digitised 1:50 000 topographic map data were provided under agreement with LINZ, Crown copyright reserved. This manuscript was reviewed internally by Simon Nathan and Peter R. Wood (IGNS) before submission, and external reviews provided by Katherine Hodgson (University of Lincoln, New Zealand) and an anonymous referee.

References Barringer, J.R.F., Lilburne, L., 1997. An evaluation of digital elevation models for upgrading New Zealand land resource inventory slope data. Proc. GeoComputation, Second Annual Conference of GeoComputation ’97 and SIRC ’97, 26^29 August 1997. University of Otago, Otago, pp. 109^116. Carrasco-Nunez, G., Vallance, J.W., Rose, W.I., 1993. A voluminous avalanche-induced lahar from Citlalte¤petl volcano, Mexico: Implications for hazard assessment. J. Volcanol. Geotherm. Res. 59, 35^46.

VOLGEO 2507 9-11-02

104

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105

Christenson, B.W., Wood, C.P., 1993. Evolution of a venthosted hydrothermal system beneath Ruapehu Crater Lake, New Zealand. Bull. Volcanol. 55, 547^565. Costa, J.E., 1997. Hydraulic modelling for lahar hazards at Cascades volcanoes. Environ. Eng. Geosci. 3, 21^30. Crandell, D.R., Miller, C.D., Glicken, H.X., Christiansen, R.L., Newhall, C.G., 1984. Catastrophic debris avalanche from ancestral Mount Shasta Volcano, California. Geology 12, 143^146. Cronin, S.J., Neall, V.E., 1997. A late Quaternary stratigraphic framework for the northeastern Ruapehu and eastern Tongariro ring plains, New Zealand. N.Z. J. Geol. Geophys. 40, 185^197. Cronin, S.J., Neall, V.E., Palmer, A.S., 1996. Geological history of the north-eastern ring plain of Ruapehu Volcano, New Zealand. Quat. Int. 34^36, 21^28. Cronin, S.J., Neall, V.E., Lecointre, J.A., Palmer, A.S., 1997a. Changes in Whangaehu river lahar characteristics during the 1995 eruption sequence, Ruapehu volcano, New Zealand. J. Volcanol. Geotherm. Res. 76, 47^61. Cronin, S.J., Hodgson, K.A., Neall, V.E., Palmer, A.S., Lecointre, J.A., 1997b. 1995 Ruapehu lahars in relation to the late Holocene lahars of Whangaehu River, New Zealand. N.Z. J. Geol. Geophys. 40, 507^520. Deutsches Zentrum fu«r Luft- und Raumfahrt, 1999. X-SAR SRTM: Shuttle Radar Topography Mission, Mapping the Earth from Space. Press and Public Relations, Cologne, 24 pp. Downs, P.W., Simon, A., 2001. Fluvial geomorphological analysis of the recruitment of large woody debris in the Yalobusha River network, Central Mississippi, USA. Geomorphology 37, 65^91. Glicken, H., 1982. Criteria for recognition of large volcanic debris avalanches. Abstr. EOS 63, 1141. Glicken, H., 1998. Rockslide-debris avalanche of May 18, 1980, Mount St. Helens Volcano, Washington. Bull. Geol. Surv. Jpn. 49, 55^106. Hackett, W.R., Houghton, B.F., 1989. A facies model for a Quaternary andesitic composite volcano, Ruapehu, New Zealand. Bull. Volcanol. 51, 51^68. Harding, D.J., Lefsky, M.A., Parker, G.G., Blair, J.B., 2001. Laser altimeter canopy height pro¢les ^ Methods and validation for closed-canopy, broadleaf forests. Remote Sens. Environ. 76, 283^297. Healy, J., Lloyd, E.F., Rushworth, D.E.H., Wood, C.P., Glover, R.B., Dibble, R.R., 1978. The eruption of Ruapehu, New Zealand on 22 June 1969. N.Z. Dept. Sci. Ind. Res. Bull. 224, 1^80. Hutchinson, M.F., 1988. Calculation of hydrologically sound digital elevation models. Third Int. Symp. on Spatial Data Handling, Sydney. International Geographical Union, Columbus, OH. Hutchinson, M.F., 1989. A new procedure for gridding elevation and stream line data with automatic removal of spurious pits. J. Hydrol. 106, 211^232. Iverson, R.M., Schilling, S.P., Vallance, J.W., 1998. Objective delineation of lahar-inundation hazard zones. Geol. Soc. Am. Bull. 110, 972^984.

Laenen, A., Hansen, R.P., 1988. Simulation of three lahars in the Mount St. Helens area, Washington, using a one-dimensional, unsteady-state stream£ow model. United States Geological Survey 88-4004. Land Information New Zealand, 1987. 1:50 000 Topographic Map, Kakatahi, Edn. 1, 1987, 260-S21. Land Information New Zealand, 1996. 1:50 000 Topographic Map, Ohakune, Edn. 2, 1996, 260-S20. Land Information New Zealand, 1996. 1:50 000 Topographic Map, Tongariro, Edn. 2, 1994, Limited Revision 1996, 260T19. Land Information New Zealand, 1997. 1:50 000 Topographic Map, Ruapehu, Edn. 2, 1994, Limited Revision 1997, 260T20. Land Information New Zealand, 1998. 1:50 000 Topographic Map, Taihape, Edn. 1, 1982, Limited Revision 1998, 260T21. Land Information New Zealand, 2000. NZ Topo Data Dictionary ^ Documentation Guide Version 2.2, unpublished. Lecointre, J.A., Neall, V.E., Palmer, A.S., 1998. Quaternary lahar stratigraphy of the western Ruapehu ring plain, New Zealand. N.Z. J. Geol. Geophys. 41, 225^246. Macedonio, G., Pareschi, M.T., 1992. Numerical simulation of some lahars from Mount St. Helens. J. Volcanol. Geotherm. Res. 54, 65^80. Madsen, S.N., Martin, J.M., Zebker, H.A., 1995. Analysis and evaluation of the NASA/JPL TOPSAR across-track interferometric SAR system. IEEE Trans. Geosci. Remote Sens. 33, 383^391. Malin, M.C., Sheridan, M.F., 1982. Computer-assisted mapping of pyroclastic surges. Science 217, 637^640. Manville, V., Hodgson, K.A., Houghton, B.F., Keys, J.R.H., White, J.D.L., 2000. Tephra, snow and water: complex sedimentary responses at an active, snow-capped stratovolcano, Ruapehu, New Zealand. Bull. Volcanol. 62, 278^293. Nairn, I.A., Wood, C.P., Hewson, C.A.Y., 1979. Phreatic eruptions of Ruapehu, April 1975. N.Z. J. Geol. Geophys. 22, 155^173. NASA/NIMA Memorandum of Understanding, summarised, 2000. Web site administrated by S. Okonek, SRTM Project O⁄ce, NASA Jet Propulsion Laboratory. http://www.jpl.nasa.gov/srtm/mou.html. Neall, V.E., Stewart, R.B., Smith, I.E.M., 1986. History and petrology of the Taranaki volcanoes. In: Smith, I.E.M. (Ed.), Late Cenozoic Volcanism in New Zealand. R. Soc. N.Z. Bull., pp. 251^263. O’Shea, B.E., 1954. Ruapehu and the Tangiwai disaster. N.Z. J. Sci. Technol. B 36, 174^189. Pairman, D., Belliss, S.E., McNeill, S.J., 1997. Terrain in£uences on SAR backscatter around Mt Taranaki, New Zealand. IEEE Trans. Geosci. Remote Sens. 35, 924^932. Palmer, B.A., 1991. Holocene lahar deposits in the Whakapapa catchment, northwestern ring plain, Ruapehu volcano (North Island, New Zealand). N.Z. J. Geol. Geophys. 34, 177^190. Palmer, B.A., Neall, V.E., 1989. The Murimotu Formation: 9500 year old deposits of a debris avalanche and associated

VOLGEO 2507 9-11-02

N.F. Stevens et al. / Journal of Volcanology and Geothermal Research 119 (2002) 89^105 lahars, Mount Ruapehu, North Island, New Zealand. N.Z. J. Geol. Geophys. 32, 477^486. Palmer, B.A., Neall, V.E., 1991. Contrasting lithofacies architecture in ring-plain deposits related to edi¢ce construction and destruction: the Quaternary Stratford and Opunake Formations, Egmont Volcano, New Zealand. Sediment. Geol. 74, 71^88. Palmer, B.A., Purves, A.M., Donoghue, S.L., 1993. Controls on accumulation of a volcaniclastic fan, Ruapehu composite volcano, New Zealand. Bull. Volcanol. 55, 176^189. Robinson, G.J., 1994. The accuracy of digital elevation models derived from digitised contour data. Photogramm. Rec. 14, 805^814. Rowland, S.K., Mackay, M.E., Garbeil, H., Mouginis-Mark, P.J., 1999. Topographic analyses of Kilauea Volcano, HawaiPi, from interferometric airborne radar. Bull. Volcanol. 61, 1^14. Schilling, S.P., 1998. LAHARZ: GIS programs for automated mapping of lahar-inundation hazard zones. US Geological Survey Open-File Report 98-638, 79 pp. Scott, K.M., Vallance, J.W., 1995. Debris £ow, debris avalanche, and £ood hazards at and downstream from Mount Rainier, Washington. US Geological Survey Professional Paper 1547, 106 pp. Shearer, J.W., 1990. The accuracy of digital terrain models. In: Petrie, G., Kennie, T.J.M. (Eds.), Terrain Modelling in Surveying and Civil Engineering. WPS, London, pp. 315^ 336. SRTM Project O⁄ce, 2001. SRTM, the mission to map the

105

world. Web site administrators E. Ramirez, S. Okonek. http://www.jpl.nasa.gov/srtm/data¢naldescriptions.html. Stevens, N.F., Wadge, G., Murray, J.B., 1999. Lava £ow volume and morphology from digitised contour maps: a case study at Mount Etna, Sicily. Geomorphology 28, 251^261. Ui, T., Kawachi, S., Neall, V.E., 1986. Fragmentation of debris avalanche material during £owage ^ evidence from the Pungarehu Formation, Mount Egmont, New Zealand. J. Volcanol. Geotherm. Res. 27, 255^264. Vallance, J.W., Scott, K.M., 1997. The Osceola Mud£ow from Mount Rainier: sedimentology and hazard implications of a huge clay-rich debris £ow. Geol. Soc. Am. Bull. 109, 143^ 163. van Zyl, J.J., 2001. The Shuttle Radar Topography Mission (SRTM): A breakthrough in remote sensing of topography. Acta Astronaut. 48, 559^565. Wahba, G., 1990. Spline models for observational data. CBMS-NSF Regional Conference Series in Applied Mathematics. Soc. Ind. Appl. Math., Philadelphia, PA. Wechsler, S.P., 2000. E¡ect of DEM Uncertainty on Topographic Parameters, DEM Scale and Terrain Evaluation. Ph.D. Thesis, State University of New York College of Environmental Science and Forestry, Syracuse, NY, 380 pp. Werner, M., 2001. Shuttle Radar Topography Mission (SRTM): mission overview. Frequenz 55, 75^79. Zebker, H.A., Madsen, S.N., Martin, J., Wheeler, K.B., Miller, T., Lou, Y.L., Alberti, G., Vetrella, S., Cucci, A., 1992. The TOPSAR interferometric radar topographic mapping instrument. IEEE Trans. Geosci. Remote Sens. 30, 933^940.

VOLGEO 2507 9-11-02