FlowDIR: a MATLAB tool for rapidly and probabilistically forecasting the travel directions of volcanic flows

Abstract

We present FlowDIR, a MATLAB tool that rapidly and objectively quantifies future travel direction probabilities for topographically controlled hazardous flows, based on analysis of summit topography. FlowDIR can achieve probabilistic forecasts of future travel directions in minutes and provides a basis for choosing the starting co-ordinates required by empirical flow models. In this work we describe the development of FlowDIR, perform a sensitivity analysis to determine the influence of input parameters on forecasted probabilities, and demonstrate its effectiveness in the retrospective forecasting of travel directions for block-and-ash flows and lava flows at three volcanoes with different summit morphologies (Shinmoedake, Colima and Merapi). In all case studies, the higher probability flow directions identified using FlowDIR agreed with the travel direction of historically observed flows. Given its intuitive outputs and rapid execution time, FlowDIR can be used to supplement existing modelling strategies for hazard assessment of topographically controlled hazardous flows prior to and during crisis. We demonstrate this by coupling FlowDIR output probabilities with an empirical hazard model to estimate probability of block-and-ash flow inundation at Gede volcano, Indonesia.

Introduction

Forecasting the areas that are likely to be impacted by hazardous phenomena is an important component of disaster risk management, and a prerequisite for pre-event land use planning and the implementation of targeted community preparedness programs (Lockwood and Hazlett 2013). Forecasting is usually achieved using numerical hazard models that emulate the natural phenomena to some extent by applying a set of assumptions and simplifications (e.g., Bonadonna 2006; Tierz et al. 2016a; Gallant et al. 2018; Clarke et al. 2020; Constantinescu et al. 2022). In this work we focus on topographically controlled hazardous flows (TCHFs), a term we use to describe volcanic flows that are generated by the extrusion and downslope movement of lava at the surface, resulting from effusive volcanic eruptions. Under this definition we consider there to be two types of TCHFs, 1) block-and-ash flows (BAFs) that are produced by the growth and collapse of a viscous lava dome or flow, and 2) lava flows, which are produced by the effusion of less viscous lavas.

A suite of tools for forecasting the spatial distribution of TCHFs currently exists; for BAFs, physics-based numerical models include Titan2D (Patra et al. 2005), VolcFlow (Kelfoun and Vargas 2016), and GPUSPH (Bilotta et al. 2016), and empirical tools include LaharZ (Schilling 1998) and the Energy Cone model (Malin and Sheridan 1982). For lava flows, physics-based models include MAGFLOW (Cappello et al. 2015) and LavaSim (Hidaka et al. 2005), and empirical models include Felpeto et al. (2001), DOWNFLOW (Favalli et al. 2005), Q-LavHa (Mossoux et al. 2016), and MrLavaLoba (de’ Michieli Vitturi and Tarquini 2018). Hazard forecasts can account for uncertainty in model inputs, and outputs (i.e., aleatoric, epistemic, or model uncertainty; Tierz 2020), by running large numbers of model simulations (typically thousands) to arrive at probabilistic forecasts (e.g. Bonadonna 2006; Tierz et al. 2016a; Gallant et al. 2018; Clarke et al. 2020; Constantinescu et al. 2022), where the frequency of the simulated scenarios represents the likelihood or the probability of the scenarios occurring (Biass et al. 2014). This can, however, be computationally demanding (Tierz et al. 2016b; Sandri et al. 2018; Clarke et al. 2020).

FlowDIR originated from the need to generate rapid probabilistic hazard footprints for TCHFs at 40 volcanoes in Southeast Asia, as part of a regional exposure assessment (see Jenkins et al. 2022). Footprint modelling requires simulation starting coordinates to be chosen, a choice that is non-trivial but that dictates which drainage(s) of a volcanic edifice TCHFs will flow down. For an individual volcano, visual examination of the volcano’s summit topography may provide insight (i.e., a flow originating from within a breached crater might be expected to travel through the breach), however, quantifying the likelihood of this scenario is prone to subjectivity, and this approach is not always feasible for a large number of volcanoes. Thus, FlowDIR was developed to forecast the expected initial travel directions of TCHFs away from a volcano’s summit, in a quantitative and objective way.

The initial travel direction, or the direction that lava extruded at a volcano’s summit moves away from the summit region (referred to as just ‘travel direction’ from hereon in), is intuitively one of the most important aspects affecting the downslope inundation area. While the transport dynamics of a lava flow depend in part on its rheological properties and effusion rate (Costa and Macedonio 2005), underlying topography is the dominant control on flow emplacement (Bilotta et al. 2019). Lava domes have been shown to exhibit preferential growth directions, which affects the direction of collapse-derived BAFs (Zorn et al. 2019). The factors that are capable of influencing dome growth directionality are not yet fully understood but mechanisms have been found to include sheer bands and slip behaviour at conduit walls (Hale and Wadge 2008), plugging of ascent pathways (Husain et al. 2014) and pre-existing topography (Walter et al. 2013).

In this work we consider topography as the major control on the travel directions of lava flows and dome collapse BAFs. We present an intuitive tool that forecasts TCHF directionality from summit topography, which can easily be applied for probabilistic forecasting. We do not attempt to supply an inundation or flow-routing model; rather, we provide a hazard assessment product that complements flow models. In the following, we describe the methodology behind FlowDIR, and retrospectively apply it to three case study examples with known flow travel directions. Following this, we assess and discuss the sensitivity of the tool, before demonstrating its application in combination with the LaharZ inundation model (Schilling 1998) for forecasting purposes. Finally, we assess and discuss the caveats.

Methods

FlowDIR is a MATLAB code that relies on the analysis of a volcano’s summit topography to estimate the likelihood of initial flow directionality for a TCHF. FlowDIR implements two complementary methods to estimate flow directionality. The first method calculates the elevation difference between the starting point and summit perimeter, it is a straight-line analysis of the topography, to quantify the relative probability of azimuthal flow from a starting point. However, since a volcano’s summit topography typically presents small-scale topographic features able to affect directionality, the second method quantifies potential flow routes using a least-cost path approach. Both of FlowDIR’s constituent methods are described in detail in Fig. 1, and in the following section, while Table 1 provides a summary of required input parameters (italicised from hereon in), along with the suggested ranges and default values.

Fig. 1

A schematic showing each step of the FlowDIR procedure: a Identification of the source; b Identification of the buffer limit; c Method 1, the azimuthal elevation difference; d Method 2, the least cost path

Table 1 A description of the FlowDIR inputs, along with a guide for their selection, suggested ranges and the default values. Default values are provided as a general guidance, however we suggest that the user becomes fully acquainted with the tools inputs before simulating

Source location

The first step within FlowDIR is shared across the two methods. The input starting coordinates are used to generate an initialisation polygon consisting of equally spaced initialisation points (Fig. 1a). Both methods are then run for each point in turn and the final result is the mean from all points. The default polygon radius is 1 DEM cell, which results in nine initialisation points. Where future vent location is not known, incorporating uncertainty into the starting location allows the user to consider multiple potential emission locations. The size of the initialisation polygon can be adjusted using the start uncertainty variable (Table 1). Although FlowDIR initialisation requires the subjective identification of a set of starting coordinates, directionality is informed based on the entire summit topography, which reduces sensitivity to the specific user-defined point. FlowDIR generates 360 radial elevation profiles at 1° azimuth increments from the starting coordinate (default length 800 m), which are referred to as swaths from hereon in (Fig. 1b). The cell with the maximum elevation along each swath is used to generate a polygon of the summit region. A buffer is then applied to this polygon to ensure that all relevant topography is included in the calculation. This is important at volcanoes with a summit dome or breached crater where the maximum elevation along the swaths oriented in the direction of the breach may be close to the initialisation point. This buffer manually extends the area included in the calculation away from the initialisation point to account for this.

Azimuthal Elevation Difference (AED)

The first method implemented in FlowDIR is called the azimuthal elevation difference (AED). Swaths are generated from each initialisation point, these consist of elevation values (in metres) that are extracted from the DEM. For each of the radial swaths (R_1:360), the elevation difference is calculated between the start and the end of the swath (Fig. 1c). Positive elevation differences show the transition from lower to higher elevations (a concave summit topography), and negative elevation differences show the transition from higher to lower elevations (a convex summit topography). Expected travel directions are characterised by relative smallest values for the elevation differences.

1º azimuth elevation differences are then averaged to wider direction bins to assess general directionality. By default, FlowDIR uses 22.5º bins with the first bin centred on north and resulting in 16 equally spaced possible directions that align with secondary intercardinal directions (i.e., NNE, NE, ENE, etc.). However, by editing the code the bin spacing can be adjusted to capture the specificity of each case-study (e.g., size of channel heads). After binning, elevation differences are inverted by subtracting each bin from the highest bin so that high elevation differences are the more likely directions. Finally, values are converted into relative probabilities by dividing each bin by the total sum of the bins (Fig. 1c). The mean of the values produced from all initialisation points within the initialisation polygon is provided for each direction bin, along with the standard error.

Although the AED approach provides a relative ranking of the likelihood of a flow to follow a specific radial direction, it conveys no information regarding the actual possibility of a flow to overtop topographic features. This is important for concave summit topographies, where very high sectors of the crater wall are unlikely to be overtopped. Therefore, as complementary information to the interpretation of radial probabilities, FlowDIR compares the binned elevation differences calculated as part of the AED procedure to a user-defined elevation threshold, and only directions with an elevation difference below the threshold are considered likely to be overtopped (Fig. 1c). This step acts as a barrier to filter out travel directions that are assigned some probability but are unrealistic due to very high crater walls. Comparison of the elevation thresholds is not included in the calculation of the AED probability, but outputs a logical response (i.e., can this sector of crater wall be overtopped yes/no). This response is depicted in the size of the diamonds for each travel direction bin in Fig. 1c, where large diamonds at a bin’s centre indicate that the crater wall in that direction is unlikely be overtopped. Setting the elevation threshold is subjective and must be explored as a function of the available knowledge about the morphology and past events at the volcano. Specifically, critical parameters to investigate include the past erupted volumes, the flow rheology, past lava heights, and previous occurrences of overtopping (Table 1). “Application to past flow travel directions” and “Application to hazard assessment: Gede volcano (Indonesia)” sections illustrate the choice of this threshold for each of our case studies.

Least Cost Path (LCP)

The least cost path (LCP) procedure maps the likely travel direction from each initialisation point using a least-cost mapping approach, where the cost raster is the inverse of the flow accumulation raster. This approach allows for material to move laterally between the DEMs elevation cells as opposed to analysing only in a specific direction along azimuth. For a given initialisation point, the underlying DEM cell becomes the first active cell. FlowDIR calculates the gradient between the active cell and the eight surrounding cells to find the one with the lowest value. This cell becomes the next active cell or the next step in the path that is taken (Fig. 1d). FlowDIR then iteratively searches for the cell with the lowest value for the gradient in the surrounding cells, making this the next active cell. At each iteration, the ID of each active cell is stored, thus adding to the path away from the initialisation point towards the edge of the summit area. The calculation is complete when the path reaches the buffer limit. To facilitate movement through the elevation matrix and to avoid the algorithm getting stuck at local minima, FlowDIR does not allow for cells to be traversed more than once for a given initialisation point: if the cell with the lowest gradient has already been active, then the 2^nd lowest cell will be used, and so on. In the case that all of the eight cells surrounding the current active cell have been traversed by the previous path of the flow, the algorithm then proceeds to the next layer of surrounding cells (i.e., the cells outside of this first eight). This is repeated until a cell that has not yet been traversed is found. For each initialisation point, FlowDIR generates a storage matrix with dimensions the same size as the DEM that records the propagation step number, where the first step closest to the initialisation point is step 1, the next step towards the buffer limit is step 2 and so on. Once the buffer limit is reached for a given initialisation point, the inverse of the step number is calculated (i.e., the intersection between the buffer limit and the path is step 1, the next cell back towards the initialisation point is step 2, and so on). Matrices for each initialisation point are stacked and summed vertically to produce a single matrix of the summed propagation step number, showing the likely flow paths from initialisation point to exit from the summit area (i.e., up to the buffer limit). Where there is more than one path out of the summit area that depends on the initialisation point used, the values of the inverse propagation step number at the buffer limit provide the likelihood of each path (Fig. 1d).

Parametric uncertainty

FlowDIR requires a minimal number of inputs (Table 1) given its basis on topography, therefore it does not need parameterisation of the flow’s physical attributes, which can be difficult to constrain prior to an eruption. However, parametric uncertainty can be explored by running a large number of simulations with parameters stochastically sampled from ranges (e.g., Tierz et al. 2016b; Gallant et al. 2018; Clarke et al. 2020; Constantinescu et al. 2022). Parameter ranges can be constrained by those provided in Table 1 and the available information about the volcano or, if running FlowDIR for multiple volcanoes, may be informed by global morphology datasets (e.g., Grosse et al. 2014). Due to its rapid execution time, exploration of the FlowDIR input parameters is achievable without the need for high performance computing facilities.

When simulating flows over natural terrain, the resolution and accuracy of the DEM is arguably one of the most important input parameters (Sparks and Aspinall 2004; Capra et al. 2011; Charbonnier et al. 2018). We suggest that the ideal DEM resolution for simulations should be a function of the topographic features that are present. Any pits or ridges in the summit area that could affect the calculation need to be resolvable. The effect of DEM resolution on outputs is tested in “Sensitivity analysis” section.

Outputs

Examples of FlowDIR outputs are provided in “Application to past flow travel directions” section, and are schematised in Fig. 1. The output from the AED procedure is a radial histogram showing the probability of travel direction (Fig. 1c) averaged over the N initialisation points, with an error bar quantifying the standard error. For each bin, diamonds show the elevation difference (in m) between the start and the edge of the crater or summit region (i.e., the raw values used to calculate the probabilities as part of the AED procedure). Diamonds are either small or large. Large diamonds indicate that the elevation difference is above the elevation threshold, which for a crater morphology suggests that flow in this direction is unlikely due to a very high crater wall, whereas smaller diamonds indicate that the elevation difference is below the threshold and flow in this direction is possible. The overtopping likelihood is only relevant to volcanoes with a concave summit, in the theoretical case of a fully convex summit, all diamonds would be the same size. The probabilities produced by the AED method can be used to populate an event tree (Newhall and Hoblitt 2002) or can be incorporated into a full probabilistic volcanic hazard assessment (PVHA) (e.g. Marzocchi et al. 2010; Sobradelo et al. 2014).

The LCP method outputs a map showing the summed inverse propagation step numbers for each cell that is traversed (Fig. 1d, “Application to past flow travel directions” section), where cells that are traversed earlier on in the path from the initialisation point to the buffer (i.e., closer to the starting coordinate) have higher step numbers than those traversed later on (i.e., closer to the buffer). This information is useful to identify the most likely route out of the summit area. To read the plot, the user should look for paths or zones of adjacent cells that extend from the starting coordinate and intersect with the buffer limit. Multiple least cost paths out of the summit area may be present due to variability between the initialisation points, in this case the most travelled route has the highest step number at the buffer limit. FlowDIR also provides the option to output iteration steps for the least cost path, and maps showing the path after every N cell steps for a given initialisation point can be produced. The outputs produced by the LCP method can complement physics- based tools to focus on the most likely sectors and to reduce the number of simulations that need running when time and/or resources are limited.

As complementary information to the results from the two procedures described above, the FlowDIR results panel also plots the DEM used for the calculation and the elevation profile of the buffer limit for reference (see “Application to past flow travel directions” section). The DEM shows obvious topographic features, such as craters or ridges, which may affect the calculation. The buffer limit profile shows any low points in the crater rim that can be used to sanity check the forecasted directionality, and to understand the effect that the crater rim has on the outputs. These plots are additionally valuable in the unlikely event that the AED and LCP methods are not in agreement, which can be an indication that the input parameters, in particular the start uncertainty and buffer width, should be re-assessed (“Sensitivity analysis” section).

Application to past flow travel directions

We applied FlowDIR to the retrospective forecasting of the directions of BAFs and lava flows at three volcanoes: Shinmoedake (Japan), Colima (Mexico) and Merapi, (Indonesia). Examples were chosen to represent a range of different summit morphologies and to include the simulation of different types of TCHFs. At the date of DEM acquisition, Shinmoedake had a wide shallow crater with a radius of ~ 400 m (20–60 m depth range). Colima had a smaller shallow crater (radius ~ 100 m, 10–27 m depth range) with no obvious breaches, whereas at the time of acquisition Merapi had a large, breached crater (radius ~ 125 m, 0–150 m depth range) (Fig. 2). For each case study we compared the directions produced using FlowDIR with the observed travel directions of subsequent historic flows that occurred after DEM acquisition, using either historical satellite imagery or Global Volcanism Program reports. Travel directions of past events are preferably described by azimuth for comparison with FlowDIR results, although many reports only provide a general direction (e.g., NW, NE). Where only a general direction is reported and satellite imagery cannot confirm a more precise azimuth, we considered the entire azimuthal bin as the directionality of the past event, e.g., for flows directed towards the south we used 169–191°, if the south bin was the highest probability produced using FlowDIR we considered this a success. In the following sections we provide more detail on the events being reproduced and the FlowDIR setup for each case study.

Fig. 2

Three-dimensional plots of the DEMs used in this study. DEM sources and resolutions are: Shinmoedake 5 m, Geospatial Information Authority of Japan; Colima 30 m, SRTM 2000; Merapi 10 m (resampled), Darmawan et al. (2017). All DEMs are plotted in the UTM coordinate system and shown at scales that best reflect the summit region

Shinmoedake (Japan)

March 2018 event description

On 1 March 2018, following several days of increased seismicity, an eruption occurred at Shinmoedake volcano within the Kirishimayama complex (Japan Meteorological Agency 2018). Japan Meteorological Agency (JMA) reported explosions on 2–3 March and on 6 March (Global Volcanism Program 2018). Explosions were accompanied by the effusion of an andesitic lava dome from the south-eastern crater (Earthquake Research Institute 2018), which gradually filled the crater floor until 9 March when lava started to overtop and flow from the NW crater wall (Earthquake Research Institute 2018). This eventually extended 200 m outside of the crater before stalling. The overtopped crater wall segment can clearly be seen in historical Google Earth satellite imagery and was measured to be 320 m in diameter at a bearing between 290–323° from the south-eastern crater (Fig. 3).

Fig. 3

FlowDIR output for forecasting the directionality of the 2018 lava flow at Shinmoedake volcano (Japan), the travel direction of the flow (the sector of the crater wall that was overtopped by lava) was between 290–323° bearing from the starting location and is shown in all panels with a red shaded area. a The DEM used for simulations (Geospatial Information Authority of Japan, 5 m resolution) with initialisation points (N = 289) shown as red crosses (please note these crosses overlap), and the crater buffer outlined in black. b The travel direction probability for each secondary intercardinal direction bin, calculated using the azimuthal elevation difference (AED) functionality. The standard error on the 289 initialisation points is shown for each bin as a red bar (very small for Shinmoedake) located at the wider end of the probability bin. Diamonds at each bin’s centre point depict the elevation difference, where larger markers show directions that exceed the elevation threshold (unlikely to be overtopped). c The least cost path (LCP) matrix, which shows the summed inverse propagation step number for each cell in the path taken from the initialisation point to the buffer limit, over all initialisation points shown in panel a, with more likely exit points indicated by higher values at the buffer limit. For Shinmoedake, due to the very flat summit topography and the high resolution DEM a highly sinuous path is taken. d The elevation profile for the buffer extent as outlined in panel a

FlowDIR parameterisation

We applied FlowDIR to assess the directionality of lava flows using a DEM from the Geospatial Information Authority of Japan (2016) (https://www.gsi.go.jp). The DEM was produced from an aerial laser survey conducted in 2016, with a spatial resolution of 5 m. In the two years between the DEM acquisition date and the 2018 eruption, another explosive eruption occurred from the south-eastern vent within the Shinmoedake crater (Global Volcanism Program 2018), but visual inspection of historical Google Earth satellite imagery suggested that this did not significantly alter the topography, meaning that the 2016 DEM was appropriate for recreating the 2018 event. We applied a buffer zone extending 50 m past the crater rim. This value was chosen based on two observations, 1) the volcano does not show a dome, and therefore does not justify using a large buffer (see “Source location” section); and 2) a small ravine visible in the DEM extended less than 50 m from the crater rim towards the SW of the crater, which may act as an outlet for a potential flow. To set the elevation threshold, we considered past events at Shinmoedake, and the topography of the crater (specifically its diameter, and the difference between the minimum and the maximum crater wall height). In 2011 a lava dome with a similar composition to the one observed in 2018 erupted at Shinmoedake (Nakada et al. 2013), forming a wide flat dome that did not sufficiently grow to flow outside of the crater. Shinmoedake’s wide, shallow crater laterally accommodates lava rather than allowing it to build tall domes. Furthermore the relatively fluid andesitic lava is unlikely to form a tall dome regardless of confinement. We selected a value of 30 m as the elevation threshold which lies in between the minimum and the maximum crater wall height. The 2018 lava dome erupted from the south-eastern vent (~ 50 × 100 m diameter) (Japan Meteorological Agency 2018), which was also the source of the explosion that occurred between DEM acquisition (in 2016) and the 2018 event. We set a starting coordinate at the centre of this south-eastern vent, in line with the starting location observed in 2018 (Earthquake Research Institute 2018). A start uncertainty of 50 m was used to reflect the diameter of the south-eastern vent, which resulted in 289 initialisation points (Fig. 3a).

FlowDIR results

Figure 3 illustrates FlowDIR outputs for Shinmoedake. Figure 3a shows the DEM; the initialisation points are marked with red crosses, and the area included in the calculation as bounded by the 50 m buffer is outlined in black. Red shaded areas show the range of bearings of the crater wall section that was overtopped by the 2018 flow measured from historical Google Earth satellite imagery (“March 2018 event description” section). Figure 3b shows the relative probability of each travel direction bin calculated using the azimuthal elevation difference (“Azimuthal Elevation Difference (AED)” section). The WNW and NW bins have the highest probabilities, with a combined relative probability of ~ 34% (Fig. 3b), which agrees with the direction of the 2018 overtopping. Error bars show the standard error computed over the 289 starting points, but the low values observed (i.e., maximum standard error of 0.05%) suggest that the AED method is not sensitive to initialisation point at Shinmoedake. Figure 3b shows that only the WNW through NW crater wall section at Shinmoedake is likely to be overtopped (small diamonds).

Figure 3c shows the results from the LCP analysis. For Shinmoedake, four points of intersection exist between the buffer limit and the summed inverse propagation step number (at 289, 300, 303, and 305°), which agrees with both the direction of overtopping in 2018 (shaded red area in Fig. 3) and the forecasted direction using the AED. In this case study, due to the flat topography of the summit crater, the path taken away from the starting coordinate (depicted in Fig. 3c by colour, where red is the start and blue is the end) is highly sinuous. The topography of the buffer extent for the starting coordinate (black line in Fig. 3a) is provided in Fig. 3d as additional information on how the crater rim elevation may affect the calculation. In the case of Shinmoedake, Fig. 3d confirms that the overtopping occurred at a topographic low in the crater wall, which is further supported by visual inspection of Fig. 2.

Colima (Mexico)

February 2002 – February 2003 event description

In 1997, Volcán de Colima entered a new phase of activity that lasted until 2011. In May 2001 a spiny lava dome began to grow within the summit crater, marking the end of more than 2 years without extrusion (Global Volcanism Program 2001). Growth continued throughout the year and, by February 2002, the dome began to overtop the crater towards the SW to produce a blocky lava flow (Global Volcanism Program 2002). From 4 February, BAFs generated by the collapse of the overtopped lava extrusion travelled towards the SW, S and W. On 9 February, BAFs travelled down the S and SW flanks of the volcano, followed on 14 February by lava flows extending from the dome up to 1.4 km towards the SW. Between February 2002—February 2003 a total of eight blocky lava flows were emplaced towards the SW and W (Global Volcanism Program 2003). Using a bin width of 22.5° this equates to a bearing of 168 – 282° from the starting coordinate (Fig. 4).

Fig. 4

FlowDIR output for forecasting the directionality of the 2002–2003 flows at Colima volcano (Mexico). Flows exited the crater in the S-W quadrant (168–282° azimuth), which is shown in all panels as the red shaded area. a The DEM (SRTM 2000, 30 m resolution) used for simulations with initialisation points (N = 25) shown as red crosses, and the crater buffer outlined in black. b The travel direction probability for each secondary intercardinal direction bin, calculated using the azimuthal elevation difference (AED) functionality. The standard error on the 25 initialisation points is shown for each bin as a red bar located at the wider end of the probability bin. Diamonds at each bin’s centre point depict the elevation difference. The presence of larger markers would indicate directions that exceed the elevation threshold (unlikely to be overtopped), however, in this case study all diamonds are the smaller size (likely to be overtopped). c The least cost path (LCP) matrix which shows the summed inverse propagation step number for each cell in the path taken from the initialisation point to the buffer limit, over all initialisation points shown in panel a, with more likely exit points indicated by higher values at the buffer limit. d The elevation profile for the buffer extent as outlined in panel a

FlowDIR parameterisation

FlowDIR was used in this case study to assess the directionality of both the lava flows and the BAFs generated between February 2002 – February 2003. We used the Satellite Radar Topography Mission (SRTM) DEM (Farr and Kobrick 2000), with an acquisition year of 2000 and a spatial resolution of 30 m. At the time of DEM acquisition Colima had only a minor crater depression (maximum depth = 27 m; Fig. 2), which was subsequently filled by the 2001 lava dome before it overwhelmed the crater resulting in lava flows and BAFs throughout 2002–2003. Using the pre-dome DEM allowed us to directly validate the ability of FlowDIR to reproduce the 2002–2003 crisis. We used the centre of the crater as the starting coordinate, in line with pictures provided in GVP reports. Given the shallow crater morphology with relatively uniform walls (i.e., 21–27 m), we did not expect any preferential crater overtopping direction. We set a elevation threshold of 30 m, greater than the maximum crater depth, to allow for all directions to be overtopped. We used a buffer of 100 m to extend the calculation outside of the summit region. Start uncertainty was set to 60 m which is the maximum distance possible that ensures all initialisation points are contained within the crater (see Table 1). Given the 30 m resolution of the DEM, this meant that two DEM cells of start uncertainty were included, reducing the over-reliance on a small number of cells.

FlowDIR results

The results of the FlowDIR analysis for Colima are shown in Fig. 4. The highest probability direction bins identified using the AED method were the SSW, SW, and WSW bins, (azimuth range 168–259°) with relative probabilities and standard errors of 10.44%/0.04% (SSW), 10.13%/0.02% (SW), and 9.50%/0.03% (WSW), respectively. All diamonds in Fig. 4c are the smaller size, meaning that the elevation change values calculated as part of the AED method are below the elevation threshold in all directions, and therefore all sections of crater wall are likely to be overtopped. This finding is in agreement with the direction of the 2002 flows (168 – 282° azimuth; marked as the red shaded area in Fig. 4). The LCP shows several paths in the SW-SE sector. The path with the greatest summed inverse propagation step numbers at the buffer limit, can be seen in green touching the buffer with an azimuth of 225–242° from the starting coordinate. This reflects the paths that are taken by the majority of the 25 initialisation points. This direction is at the centre of the direction range taken by the 2002–2003 flows and agrees with the highest probability travel direction bin (SW) calculated with the AED method. FlowDIR performed well for this case study since both methods were able to identify the predominant travel direction.

Merapi (Indonesia)

December 2018 – September 2019 event description

In May–June 2018, after several years of quiescence, phreatic explosions occurred at Merapi, followed on 12 August 2018 by the effusion of a lava dome at the summit (Global Volcanism Program 2019). Dome growth direction was towards the SE crater wall breach (Kelfoun et al. 2021) generated during the 2006 eruption (Ratdomopurbo et al. 2013). By December, ‘block avalanches’ (based on the context within reports, these are assumed to be sizeable rockfalls and not BAFs since those are described separately) were reported up to 1 km away in the Gendol drainage on the SSE flank of the volcano which is directly fed by the SE breach in the crater (oriented between 100–205° from the crater centre as measured in historical Google Earth satellite imagery). These block avalanches continued throughout the start of 2019 (Global Volcanism Program 2019). On 18 February 2019 BAFs travelled up to 1 km down the Gendol drainage, continuing into September accompanying low volume dome growth. For forecasting this event we consider our simulations a success if the highest probability travel direction lies within the 100–205° oriented breach (Fig. 5).

Fig. 5

FlowDIR output for forecasting the directionality of the 2018–2019 BAFs at Merapi volcano (Indonesia), the direction of the flow is shown in all panels with a red shaded area. a The DEM (Darmawan et al. 2017, 10 m resolution) with initialisation points (N = 49) shown as red crosses, and the crater buffer outlined in black. b The travel direction probability for each secondary intercardinal direction bin, calculated using the azimuthal elevation difference (AED) functionality. The standard error on the 49 initialisation points is shown for each bin as a red bar located at the wider end of the probability bin. Diamonds at each bin’s centre point depict the elevation difference. The presence of larger markers would indicate directions that exceed the elevation threshold (unlikely to be overtopped), however, in this case study all diamonds are the smaller size (likely to be overtopped). c The least cost path (LCP) matrix which shows the summed inverse propagation step number for each cell in the path taken from the initialisation point to the buffer limit, over all initialisation points shown in panel a, with more likely exit points indicated by higher values at the buffer limit. d The elevation profile for the buffer extent as outlined in panel a

FlowDIR parameterisation

We used the DEM of Darmawan et al. (2017), generated from a combination of terrestrial laser scanning (TLS) and structure-from-motion from unmanned aerial vehicles (UAV). The TLS campaign was conducted in September 2014, while the UAV flights were in October 2015 meaning that the DEM acquisition occurred during the period of quiescence that preceded the appearance of the lava dome in 2018. The initial 0.5 m resolution DEM was re-interpolated to 10 m using a Nearest Neighbour method to speed up computations. At DEM acquisition date, Merapi had a large, breached crater that was open towards the SSE. Photos included in the Global Volcanism Program bulletins show that the 2018 lava dome started to grow from a central location within the crater, which we used here as the starting coordinate. The DEM also shows a variable crater wall height above the base, ranging from 0 m in the SE where the crater is breached, to ~ 70 m in the NW and ~ 150 m in the NE (Fig. 5d). Given this disparity in the heights of the different sectors we expected that overtopping of the higher sectors is unlikely. We therefore set the elevation threshold to 100 m, a value that limits overtopping in the NE sector. A start uncertainty of 30 m was chosen to cover the crater floor. Merapi has a breached crater with an older dome at the summit, this means that for swaths in the direction of the breach the highest point along the swath is the starting point. Therefore, to ensure inclusion of the dome topography in the direction of the breach, we applied a buffer of 100 m, which is the approximate distance from the starting coordinate to the crater wall.

FlowDIR results

Figure 5b shows that most likely directions were, in decreasing order, SSE, S, and SE with relative probabilities and standard errors of 11.76%/0.08% (SSE), 10.39%/0.2% (S), and 10.17%/0.09% (SE) respectively. The low standard error (shown as red lines on each bin) reveals that there was little variability between the 49 initialisation points. The elevation change in each direction was less than the 100 m threshold we set (as with Fig. 4b, all diamonds in Fig. 5b are small), suggesting that all directions could be overtopped. Figure 5c shows the least cost paths taken from the start to the buffer. There are several paths, three of which are within the SE breach (126, 144, 172°), which is in agreement with the direction taken during the 2018–2019 eruption and one which extends out of the crater in the WNW direction (282°). The colour map of the pathways indicates that the WNW oriented path (green colour) has the higher summed inverse propagation step number at the intersection with the buffer limit, and was therefore the more traversed. Figure 2 shows that the WNW direction has a considerable secondary breach in the crater wall, therefore even though this travel direction was not taken during the actual eruption, there was a high chance that it could have been according to our analysis.

Sensitivity analysis

We performed a sensitivity analysis of each of FlowDIR’s three main input parameters (the buffer width, the start uncertainty, and the DEM resolution) on the estimation of flow directionality using both methods (AED and LCP) with a one-at-a-time approach (i.e., varying one parameter whilst keeping the others fixed to the values in Table 2). We did not consider the elevation threshold since it generates a binary response. For the buffer parameter, which should be kept as low as possible while including all relevant topography outside the crater (e.g., breaches or smaller channels in the crater wall; “Source location” section), we tested a range between 20–200 m (Table 3). For the start uncertainty, we tested a range between 10–50 m (Table 3), meaning that initialisation points were simulated for crater widths up to 100 m. This upper bound was chosen to limit the run times since 50 m of start uncertainty resulted in 289 initialisation points for the Shinmoedake 5 m DEM. We varied the DEM resolution between 5–30 m (Table 3) (using Nearest Neighbour interpolation), ranging from the maximum resolution advisable for use in FlowDIR regarding computation time (“Computation time” section) to the minimum resolution advised to capture topographic features relevant to flow direction. Since upsampling of the DEM does not result in the addition of new information, Colima was not considered.

Table 2 FlowDIR inputs used for recreating past events at Shinmoedake (March 2018), Colima (February 2002—February 2003), and Merapi (December 2018 – September 2019) volcanoes

Table 3 Parameters for the sensitivity analysis of FlowDIR applied to the three case studies. When not being tested, parameters were set to the values used in the case studies (Table 2). Note that due to differences in the DEM resolutions between volcanoes the start uncertainty measured in metres equates to a different number of DEM cells for the case studies

To test the sensitivity of the AED method output variables to the relevant inputs, we evaluated change in the distribution of output probabilities. For ease of interpretability we examined eight direction bins as opposed to the 16 calculated and shown in Figs. 3, 4 and 5. For the LCP method, we measured the azimuths of the points of intersection between the least cost paths and the buffer limit.

Azimuthal elevation difference

Figures 6a-c show the sensitivity of the radial probability for 8 quadrants as a function of the buffer width, start uncertainty, and DEM resolution where flat lines indicate low sensitivity. The AED method is moderately sensitive to the width of the buffer for Shinmoedake, Colima and Merapi (i.e., maximum change of 4.75% for the W bin at Colima volcano when the buffer is increased from 100 to 200 m) (Fig. 6a). This is expected given that the AED calculation is based on the elevation difference between the initialisation point and buffer limit. This indicates that the buffer width should not extend to unnecessary distances outside of the summit region at the risk of incorporating irrelevant topography into the calculation. Figure 6b shows the sensitivity of FlowDIR’s AED method to the start uncertainty parameter. All case studies show low sensitivity to this parameter, with consistent probabilities of flow direction across the ranges of start uncertainty tested. DEM resolution was found to have a minimal impact on the AED (Fig. 6c).

Fig. 6

Plots showing the sensitivity of FlowDIR outputs to buffer length (a, d), uncertainty on starting coordinate (b, e) and DEM resolution (c, f) for a-c) the AED method and d-f) the LCP method. The AED method illustrates the relative probability of each direction bin with variable input parameters. For ease of interpretability we examined eight direction bins as opposed to the 16 calculated. The LCP method illustrates the azimuth of paths that intersect the buffer limit (i.e., paths from the starting coordinate out of the crater). The grey dashed lines in figures d-f show the limits of the 16 bins. Steeper lines at varying values indicate higher sensitivity

Least cost path

Figure 6d-f show the sensitivity of the LCP method to the input parameters. Markers that are horizontally aligned across the range of inputs tested suggest low sensitivity. In some cases, a change in the input parameter results in the addition or removal of new paths out of the summit area (i.e., new points of intersection with the buffer). Figure 6d shows that the LCP method has limited sensitivity to the buffer width across the range tested. For Shinmoedake, all tested buffer values consistently result in four points of intersection with the buffer that lie between 286–307°. For Colima, the majority of points show a spread between 124–236° but the northern paths are identified at a lower range of buffer values (Fig. 6d). This larger range is expected given the coarse resolution of the DEM. Nevertheless, this range still results in the majority of paths to be identified within the same sector of the volcano, which we consider as an acceptable level of sensitivity. For Merapi, two directions are taken across the range of buffers tested (see Fig. 5c for reference): one path towards the WNW (279–289°) and a set of paths out of the breach (118–199°) (Fig. 6d). The wide azimuth range taken by paths in the direction of the breach reflects the size of the breach, illustrating a variability associated with the choice the initialisation point for individual simulations.

The sensitivity of the LCP to the start uncertainty input parameter is negligible for Shinmoedake (Fig. 6e). For Colima, at the lower values of start uncertainty tested (10–40 m), only one path located at 233° is taken towards the buffer, whereas a start uncertainty of 50 m results in a spread in output azimuths between 134–233°. For Merapi, two routes out of the summit area (WNW, 282–283°, and the SE breach, 125–173°) are consistently displayed. However, the SE breach oriented paths show a wider spread with increasing start uncertainty.

The LCP applied to Shinmoedake shows sensitivity to the DEM resolution, with a coarse resolution resulting in a wider range of identified paths (Fig. 6f). For a resolution of 10 m paths also occur NE of the starting coordinates (at 51 and 57°). For Merapi, reducing the DEM resolution from 10 to 30 m results in a loss of the identification of the WNW oriented path (282°). The sensitivity of results to the DEM resolution is a function of the topography itself and is expected to differ between case-studies. Nevertheless, we suggest that an ideal resolution of 5–10 m should be used in order to resolve relevant topographic features for TCHF.

Computation time

All simulations were run using a MacBook Pro 2021 laptop (Apple M1 pro chip with a 10-core processor and 16 GB RAM). We found that the computation time needed to run one full simulation (which consists of both the AED and the LCP) was dependent on both the input parameters and on the topography of the case study (Fig. 7). The width of the buffer had no impact on the simulation time across the range tested (black lines in Fig. 7), while increasing the start uncertainty resulted in longer simulation times (blue lines in Fig. 7). Coarser DEM resolution resulted in shorter simulation times for Shinmoedake (red lines in Fig. 7). Using Shinmoedake as an example, the DEM resolution parameter had the most influence on the computation time (0.2 min for 30 m and 60 min for 5 m) while for the start uncertainty parameter the computation time varied between 5–60 min. Simulation times for Shinmoedake far exceeded the time required to run simulations for the other volcanoes, which is in part due to the high DEM resolution (5 m), the choice of baseline parameters (higher end of the start uncertainty range), and the topography. The LCP method will continue simulating in the direction of least cost until the buffer is reached. Figure 3c illustrates that the wide flat crater (~ 400 m radius) necessitated almost 1,000 steps to reach the buffer, contributing towards the prolonged simulation times. While drone photogrammetry has made high resolution DEMs increasingly accessible, the increased simulation time required to run FlowDIR using a fine resolution DEM – with marginal to no improvement in output – means that DEMs below 5 m resolution should be reserved for situations where the problem scale necessitates it, such as small-scale volcanic features like cinder cones. Volcanoes have a range of summit topographies, and simulation times are more likely to be towards the lower end of the range demonstrated here. For a more rapid result, uncertainty in the start point might be reduced, however this should be done in the context of the topography, and the use of some start uncertainty is recommended.

Fig. 7

Computation time of FlowDIR as a function of: the buffer width (solid black lines), DEM resolution (dotted red lines) and uncertainty (dashed blue lines). When not the parameter being tested, the values of the remaining parameters were set to those used in the case studies and detailed in Table 2. The original resolution of DEMs used are: Shinmoedake 5 m; Colima 30 m, Merapi 10 m

Application to hazard assessment: Gede volcano (Indonesia)

FlowDIR can be used to complement existing strategies for TCHF hazard assessment. By coupling FlowDIR with a simple empirical model (e.g., Energy Cone or LaharZ), uncertainty can be incorporated into an otherwise deterministic approach with little further computational expense. We illustrated this here by applying FlowDIR to the assessment of BAF hazards at Gede volcano (Indonesia) by combining the output of the AED method with a BAF calibrated version of LaharZ (after Widiwijayanti et al. 2009). We used the Indonesian national DEM (DEMNAS: Geospatial Information Agency, 2018) acquired between 2010–2015 at 8.3 m resolution. The summit of Gede consists of a 900 m wide crater area (0–170 m depth range) that is filled with an andesitic lava dome erupted in 1840 (Belousov et al. 2015). We set the FlowDIR starting coordinate to the centre of the eroded 1840 dome under the assumption that the location of future dome growth will follow past events. We used a start uncertainty of 50 m (i.e., 169 individual points) to ensure that the majority of the dome was covered by initialisation points (Fig. 8a). Gede’s wide crater is breached towards the north, meaning that the highest elevation along all swaths in this direction is at the start of the swath. Therefore, to ensure that the calculation extended away from the starting coordinate in the direction of this breach we used a buffer of 100 m. This was a trade-off between ensuring that the northern topography was analysed, while limiting the amount of the southern flank past the crater rim to be included. As the southern crater wall at Gede rises 170 m above the crater floor, and given the disparity between the height of the northern (breached) and southern sectors, we set the elevation threshold to 100 m, meaning that only the highest section of crater wall cannot be overtopped.

Fig. 8

Coupling of FlowDIR output coordinates and probabilities produced using the azimuthal elevation difference to initialise a BAF calibrated version of the LaharZ code with an input volume of ${9.8 \times 10}^{6} {\text{m}}^{3}$ at Gede volcano (Indonesia). a Zoomed view of the summit showing the buffer limit, initialisation points and the LaharZ starting coordinates. b Probability of inundation produced by running simulations from the coordinates shown in a. A buffer of 300 m was added to LaharZ footprints as the minimum expected channel overtopping due to the surge component (after Widiwijayanti et al. 2009). c Azimuthal elevation difference output from FlowDIR. Coordinate reference system: WGS84/UTM zone 48S

The AED method outputs the UTM coordinates of each bin’s central swath at the buffer limit, which we used as LaharZ starting coordinates, omitting bins where overtopping was unlikely (Fig. 8a,c). For each set of LaharZ starting coordinates we simulated a flow volume of ${9.8 \times 10}^{6} {\text{m}}^{3}$, representing the ${90}^{\text{th}}$ percentile of BAF volumes in the FlowDat global dataset, (Ogburn 2016). Each LaharZ output footprint was extended laterally by 300 m to represent the possibility of the surge component escaping the channel (Widiwijayanti et al. 2009; Lerner et al. 2022). Each flow inundation footprint was then weighted by the respective travel direction probabilities in its respective bin (Fig. 8c). Finally, rasters were aggregated to produce a probabilistic BAF inundation map for Gede (Fig. 8b).

We found that, for Gede volcano, the NE and SW flanks of the volcano were the most likely to be affected by BAFs (conditional probabilities of inundation were 0.68 and 0.16 respectively, Fig. 8b). This agrees with the directions of past events at Gede (Belousov et al. 2015; Tennant et al. 2021).

Outputs such as Fig. 8 provide relative probabilities of hazard inundation, which are conditional to the occurrence of the flow. This may be incorporated into a PVHA framework that not only examines the probability of inundation but looks at the absolute probability of the event and its preceding events happening in a given time window. These frameworks typically take the form of an event tree, which consists of increasingly specific nodes or events branching from a necessary prior event. Nodes are often represented as probability distributions, and in a Bayesian framework, prior information about a system can be updated as new information becomes available. Tools are available to assist with Bayesian PVHA (BET_VH, Marzocchi et al. 2010; HASSET, Sobradelo et al. 2014). These tools require the user to input hazard information, which in BET_VH takes the format of a hazard raster, while in HASSET this is a probability of inundation per sector of the volcano. FlowDIR can assist with this; as a stand-alone tool it may be used to provide probabilities of inundation per sector, or, as demonstrated in Fig. 8 when coupled with an inundation model, rasters may be generated for incorporation into BET_VH.

Considerations

The Shinmoedake case study showed FlowDIR to be highly accurate in forecasting the azimuth while using a high-resolution DEM (5m). While the Colima example demonstrated FlowDIR’s ability to forecast the most likely travel direction when this was not apparent through visual examination of the DEM (i.e., all crater walls sections are of similar height). The Merapi case study highlighted the value of quantitative approaches to forecasting travel directions since both the SSE and WNW directions had high FlowDIR forecasted probabilities suggesting that, while not actualised in the 2018–2019 eruption, flows directed towards the WNW were likely. The sensitivity of FlowDIR results to the choice of input parameters varied between the case studies, which reflected the variability in the topography between the volcanoes and the differences in the DEM resolution. Moderate sensitivity of outputs to the buffer width parameter highlights the importance of including this in simulations such that excessive topography is not included in the calculation. This also emphasizes the value of considering all FlowDIR outputs together.

While FlowDIR has been proven to capture TCHF flow direction well, the following caveats should be considered before use:

• FlowDIR only considers topography to estimate travel direction. This approach is commonly used in TCHF simulations (e.g., Felpeto et al. 2001; Favalli et al. 2005; de’ Michieli Vitturi and Tarquini 2018), and allows computation time to be significantly reduced, enabling probabilistic analyses.

• Simulations are heavily reliant on the availability of an up-to-date DEM.However, active craters continuously undergo a combination of erosive and depositional processes that dynamically alter topography. FlowDIR therefore captures the probability of travel directionality at a fixed point in time, and any change to topography requires FlowDIR outputs to be carefully and critically discussed. With the rising use of unmanned aerial vehicles for monitoring volcanoes (James et al. 2020) the generation of new DEMs during or immediately after eruptions is becoming increasingly common using LiDAR or photogrammetry techniques (e.g., De Beni et al. 2019; Biass et al. 2019; Walter et al. 2020; Román et al. 2022; Civico et al. 2022).

• FlowDIR is suitable for forecasting the directionality of BAFs and lava flows. However, for lava flows this is limited to flows extruded from a central vent, so the tool is not suitable yet for forecasting the directionality of flows originating from fissure eruptions. In FlowDIR the source term is defined by a square, whereas for a fissure eruption the source can be highly elongated (Kīlauea’s 2018 lower East Rift Zone fissure extended 6.8 km in length; Neal et al. 2018; Meredith et al. 2022).

• BAFs can be produced by a variety of different dome collapse mechanisms including gravitational loading (Ui et al. 1999), internal gas overpressures (Voight and Elsworth., 2000), topography-controlled collapses (Voight et al. 2002), intense rainfall events (Carn et al. 2004), and changing extrusion directions (Loughlin et al. 2010). Given its basis on topography, FlowDIR should be limited to forecasting the directionality of BAFs produced by topography-controlled collapses (which occur when the dome exceeds the crater walls), and gravitational loading (which includes oversteepening) (Harnett et al. 2019). While it is not possible to forecast the specific mechanism that will lead to dome-collapse at a volcano, topography and gravitational loading have accounted for 66% of the globally catalogued dome collapse events (Harnett et al. 2019), meaning that FlowDIR has a wide scope for forecasting applications.

• As part of the azimuthal elevation difference method, swaths are binned into azimuthal directions relative to the starting point. The assignment of bins can be influenced by the choice of starting points. For instance, cells initially categorised in the North bin may end up in the South bin if the starting point is shifted Northward. However, this issue is typically relevant only in simulations conducted at volcanoes with wide craters, a high-resolution DEM, and a large amount of start uncertainty.

Conclusions

Knowledge about the expected travel directions of topographically controlled hazardous flows is essential information required for the reduction of risk to communities living around active volcanoes. Tools used for probabilistic hazard forecasting can be challenging to set-up and use, and can be computationally, data, and time intensive. To simplify matters, we have developed FlowDIR, an open-source, user-friendly topography-based approach for the rapid probabilistic assessment of initial travel directions. Applications of FlowDIR to the retrospective forecasting of past travel directions have demonstrated its ability to reproduce the past events at three volcanoes (Shinmoedake, Colima and Merapi) with differing summit topographies and digital elevation basemaps. FlowDIR results can be combined with empirical hazard models to produce probabilistic estimates of TCHF inundation areas with little computational expense. For more computationally demanding physics-based hazard modelling of flows, FlowDIR can be used to prioritise simulations and provide an evidence base for identifying more likely paths taken away from the summit area. Alternatively, if a volcano starts to show signs of unrest, FlowDIR can be used as a stand-alone tool to provide a rapid forecast of the flank that is likely to be impacted.

To date, FlowDIR has been tested on four volcanoes with various summit topographies. Future work would benefit from further validation with additional case studies encouraging collaborative development within the hazard community. With this in mind, FlowDIR has the potential to become a practical and useful tool for hazard modellers, government agencies, researchers, and volcano observatory staff alike.