The influence of surficial features in lava flow modelling

Lava flows can cause substantial and immediate damage to the built environment and affect the economy and society over days through to decades. Lava flow modelling can be undertaken to help stakeholders prepare for and respond to lava flow crises. Traditionally, lava flow modelling is conducted on a digital elevation model, but this type of representation of the surface may not be appropriate for all settings. Indeed, we suggest that in urban areas a digital surface model may more accurately capture all of the obstacles a lava flow would encounter. We use three effusive eruption scenarios in the well-studied Auckland Volcanic Field (New Zealand) to demonstrate the difference between modelling on an elevation model versus on a surface model. The influence of surficial features on lava flow modelling results is quantified using a modified Jaccard coefficient. For the scenario in the most urbanised environment, the Jaccard coefficient is 40%, indicating less than half of the footprints overlap, while for the scenario in the least urbanised environment, the Jaccard coefficient is 90%, indicating substantial overlap. We find that manmade surficial features can influence the hazard posed by lava flows and that a digital surface model may be more applicable in highly modified environments.


Introduction
Recent effusive volcanic eruptions such as the 2018 Kīlauea eruption in Hawaii, USA, have reminded the global community how disruptive lava flows can be. Although few deaths are attributed to lava flows (Harris 2015;Brown et al. 2017), they can have severe, immediate impacts to the built environment and prolonged economic and societal consequences. For example, the 1973 Vestmannaeyjar Volcanic Field eruption on the Icelandic island of Heimaey threatened the local fishing harbour (e.g. Williams and Moore 1983;McPhee 1989) Bonaccorsco et al. 2016;Rongo et al. 2016). In yet another example, over the course of 3 months in 2018, the development of the lava flow field on Kīlauea's Lower East Rift Zone (Hawaii, USA) destroyed over 700 houses (Neal et al. 2019) and the local government has reported small business closures over a year later (County of Hawai'i 2019). These consequences have been recognised for decades. In fact, Thomas Jaggar devised ways to protect Hilo Harbour (Hawaii, USA) in future eruptions in the 1940s (Jaggar 1945). Such planning continues to this day (CDEM 2015). As with all hazards, precrisis planning can minimise the psychological and physical impacts and recovery costs for communities affected by lava flows (UNDDR 2019). Preparation and mitigation actions can take many forms, such as the development of plans and policies (e.g. the creation of the Auckland Volcanic Field contingency and evacuation plans (CDEM 2015;Wild et al. 2019)), lava flow modelling (e.g. the DOWNFLOW modelling undertaken by Favalli et al. (2006Favalli et al. ( , 2009a and Chirico et al. (2009) in Goma, Democratic Republic of Congo), volunteer and professional trainings (e.g. the electricity utility on Hawai'i Island creating trainings about how to strengthen their networks to withstand the high temperatures of lava flows ), and exercises (e.g. Exercise Ruamoko, a New Zealand all-of-nation desktop exercise focused on unrest in the Auckland Volcanic Field (Brunsdon and Park 2009;Lindsay et al. 2010)). All these activities rely on understanding the potential hazard(s), in this case when and where lava flows can occur.
Assessing the potential lava flow hazard to a community can be difficult for a variety of reasons. Many volcanoes have not been studied in detail, and even well-studied volcanoes commonly lack extensive research on specific characteristics of previous individual eruptions; metrics such as recurrence interval or average lava flow episode volume are often unknown (e.g. Mt. Cameroon, Cameroon; Bonne et al. 2008;Favalli et al. 2012). Additionally, predicting the location of the next volcanic vent is critical to predicting lava flow inundation areas (Connor et al. 2012;Bilotta et al. 2019); this can be especially challenging in volcanic fields where the area under consideration for vent opening may be particularly large with few spatio-temporal trends (e.g. Allen and Smith 1994;Gallant et al. 2018).
In some volcanic regions the natural topography has been anthropologically altered during urbanisation (e.g. Al-Madinah, Saudi Arabia (Runge 2015), Auckland, New Zealand (von Hochstetter 1859;Searle 1964)). In these cases, the increasing number of obstacles (i.e. elements in the built environment) a lava flow may encounter may be increased; in other words, lava flows will not simply follow the natural topography . Despite substantial anthropogenic alteration to some natural topographies, lava flow modelling is still normally conducted on digital elevation models (DEMs) that do not consider such alterations (e.g. Kereszturi et al. 2014).
Digital elevation models are commonly used when modelling volcanic hazards. One modern method to create a DEM is using a LAS dataset, which is created using a Light (intensity) Detection and Ranging (LiDAR) system (e.g. White and Wang 2003;Favalli et al. 2009b). To create the dataset, a LiDAR system is flown over the area of interest while a laser in ultraviolet, visible, or infrared wavelengths is directed at the ground (e.g. Wagner et al. 2006;Heidemann 2018). As the LiDAR system flies over the area of interest, it measures the reflectance of the original wave emitted (e.g. Axelsson 1999;Wagner et al. 2006;Favalli et al. 2009b). Several returns are measured (e.g. Wagner et al. 2006), with each return providing different information about the ground below. The returns are numbered based on the order in which they are sensed (e.g. Heidemann 2018). The number of returns measured depends on the system being used, but the early returns detail surficial features while later returns describe lower features (e.g. the bare ground; Axelsson 1999). The first return is commonly used to create a digital surface model (DSM), which includes the bare ground and surficial features including trees, buildings, etc. (e.g. Wagner et al. 2006;Heidemann 2018). One of the later returns is used to create a DEM. In-between returns can be used to create hybrid surface models (e.g. a DSM that does not include trees) through the process of classification or filtration (e.g. Axelsson 1999;Wagner et al. 2006). A second modern method to create a DEM is using Structure-from-Motion photogrammetry (e.g. Westoby et al. 2012;Ouedraogo et al. 2014). In this method, a set of overlapping but offset photographs is taken (e.g. Westoby et al. 2012;Dietrich 2015). Algorithms can then identify features from multiple images at different angles to create a three-dimensional surface. Sufficient images are taken to blanket an area with images, enabling a DEM (or DSM) to be created (e.g. Westoby et al. 2012;Ouedraogo et al. 2014;Dietrich 2015).
Here, we present lava flow hazard modelling for three effusive eruption scenarios in the Auckland Volcanic Field (AVF), which underlies the city of Auckland, New Zealand, to highlight the effects of using different elevation models, namely DEMs and DSMs. In Section 2, we provide an overview of the AVF and the scenarios used, which were created in the context of the Determining Volcanic Risk in Auckland (DEVORA) research programme. In Section 3, we describe our methods, including how we selected the MOLASSES lava flow model and how we compared our modelling results. In Section 4, we present the MOLASSES results from modelling the DEVORA Scenarios on both a DEM and a DSM. Finally, in Section 5, we detail some of the implications of modelling lava flows on a DEM versus a DSM and compare our modelling results to the lava flow hazard footprints developed in an earlier version of the DEVORA scenarios.

Auckland volcanic field, North Island, New Zealand
The AVF is a monogenetic, basaltic volcanic field situated on the North Island of New Zealand (Fig. 1). It has been active for approximately 190kyr with the most recent eruption 550 yr BP (Leonard et al. 2017). Since the last eruption, New Zealand's largest city, Auckland, has been built on top of the volcanic field, and consequently, Auckland could be severely impacted by a future AVF eruption (Deligne et al. 2015;Hayes et al. 2017). There are neither strong spatial-temporal nor spatialvolumetric trends in the AVF (Bebbington and Cronin 2011;Bebbington 2015), although there are trends correlating bulk rock chemistry with distance to neighbouring vents, volume, and time (Le Corvec et al. 2013;McGee et al. 2015). Without clear spatial-temporal or spatial-volumetric trends (Bebbington and Cronin 2011;Bebbington 2015), it is difficult to assess probabilistic hazard for the city, including location and probability of the next vent opening (Lindsay et al. 2010;Sandri et al. 2012;Hayes et al. 2018). The Determining Volcanic Risk in Auckland (DEVORA) research programme (http:// www.devora.org.nz/) was established in 2008 to study the volcanic risk in the city from local and distal volcanoes. To aid stakeholders in preparing for an eruption, the DEVORA research team developed eight hypothetical AVF eruption scenarios (Hayes et al. 2018(Hayes et al. , 2020, four of which include lava flows ( Fig. 1).

DEVORA scenarios
The DEVORA scenarios are hypothetical eruption sequences, incorporating multi-hazard modelling (Hayes et al. 2018(Hayes et al. , 2020. The hazards were modelled to represent the full range of possible eruptive phenomena and hazard intensities (Hayes et al. 2018(Hayes et al. , 2020. Not all hazards are included in every scenario; rather, hazards are included based on how frequently they have occurred in past AVF eruptions, and the local environmental conditions at the scenario vent area (e.g. onshore/offshore; depth to water table). The resulting DEVORA scenarios have been or will be used in a variety of projects conducted by DEVORA researchers and their stakeholders, including calculating the impact of an AVF eruption on electricity transmission networks (Tsang SWR, Lindsay JM, Kennedy B, Deligne: Thermal impacts of basaltic lava flows to buried infrastructure: workflow to determine the hazard, submitted), reviewing the local government volcanic contingency plan (A. Doherty, pers. comm.), assessing waste disposal during the recovery phase of an AVF eruption ( In these scenarios, lava flow hazard was represented by a temporal series of hand-drawn areal lava flow footprints, following the precedent of hand-drawing an areal footprint set by Deligne et al. (2015). In Hayes et al. (2018), topography maps were used to estimate the steepest gradient at the flow front's location. The flow front was then drawn to advance a short direction in the direction of the steepest gradient. Once the footprint had been drawn, the volume represented on the map was calculated using the area of the flow's footprint and an average (i.e. mean) thickness based on analogue lava flows. This process was used until the total volume prescribed by the scenario had been distributed. One of the motivations of our study was to update the hand-drawn lava flow hazard footprints in these scenarios using numerical modelling.

Methods
To analyse the influence of the built environment on lava flow modelling, we quantitatively model the lava flow hazard for three of the DEVORA Scenarios (Mt Eden, Birkenhead, and Ōtāhuhu; Table 1) on both a DEM and a DSM. Lava flow modelling of the Rangitoto Island scenario was not undertaken because the scenario vent is within 25 m of the Hauraki Gulf, minimising the possibility of built environment impacts. Although the vent in the Ōtāhuhu scenario also lies at the edge of a tidal area, the closest body of water is very shallow.

Lava flow hazard model selection
As our modelling results are intended to replace the hand-drawn hazard footprints presented in Hayes et al. (2018), our model selection was guided by the lava flow characteristics included in the original scenarios together with DEVORA stakeholder and research programme guidelines in place for selecting volcanic hazard models for the Auckland Volcanic Field.
The selection of the most appropriate lava flow model(s) is vital to ensure the necessary outputs for the intended application. Hayes et al. (2018) includes quantitative descriptions of the lava flow's presence (areal footprint) and thickness. Seven published and peer-reviewed models produce both of these outputs: LAVASIM, (See figure on previous page.) Fig. 1 Map showing Auckland with the volcanic centres of the AVF. The black oval shows the possible boundary of the AVF based on the known eruptive centres, shown as black triangles (Runge 2015;Kereszturi et al. 2014). The DEVORA scenario vents that we use as case studies are denoted with solid, red triangles; the letters in each triangle indicate the scenario. B is for Birkenhead; the Birkenhead scenario has two vents, both of which fall within the northernmost red triangle. M is for Mount Eden; O is for Ōtāhuhu. Inset shows the North Island of New Zealand and major cities; Auckland is boxed. Figure modified from Bebbington (2015) and Hayes et al. (2018). Data from Kermode (1992) and Land Information New Zealand MAGFLOW, MOLASSES, MULTIFLOW, SCIARA, and VOLCFLOW. Additionally, the DEVORA research team use the following criteria to select hazard models used for the AVF: Ideally hazard models should be explored that can be run in-house in a reasonably short timeframe should an eruption begin. Codes should ideally be open access, i.e. not dependent on a paid software license. If necessary, models that are dependent on proprietary software (e.g. ArcGIS) can be considered if the software is widely available to researchers and stakeholders. Since we are modelling processes in a highly topographically modified and densely populated area, codes that can handle data with high spatial resolutions are preferred. Preference is given to codes that do not need to be significantly modified to be applied to Auckland context. The code(s) (with minimal or no modification) must be executable on a desktop computer. Finally, there is a preference for codes that have already been validated on natural case studies.
Based on these guidelines, we selected the MOLASSES lava flow simulation code (Gallant et al. 2018). MOLAS-SES determines an areal footprint and the thickness of the lava flow at given locations, yielding the required outputs for our case studies, and can be run locally on both a DSM and DEM without requiring a proprietary license. For a more detailed explanation of how and why MOLASSES was selected, see Tsang (2020).

Lava flow hazard modelling
We have run MOLASSES six times here to illustrate how results can differ if a surface model is used instead of an elevation model. All three scenarios were run twice in MOLASSES, once on a DEM and once on a DSM. The DEM and DSM used in this project are publicly available from Land Information New Zealand (LINZ) and were created using the same LAS dataset captured by Auckland Council during the second half of 2013. They both have a resolution of 1 m and a vertical accuracy of ±0.2 m and a horizontal accuracy of ±0.6 m (both at a 95% confidence level). The same source DEM and DSM were used for all model runs although the DEM and DSM were cropped to the area of interest in each scenario. No special processing was applied to either product; thus, the DSM includes all surficial features, including trees. In our DSM, many of the buildings in the localised areas of interest are used as housing and generally vary from approximately the lava flow's mean thickness to many metres above the lava flow's upper crust. The input parameters for the MOLASSES modelling are included in Table 2. Since the lava flow modelling results presented in this paper effectively update the qualitative footprints in Hayes et al. (2018), the lava flow hazard outputs presented should be referred to as the DEVORA scenarios, version 1.1.

Quantifying the influence of surficial features on lava flow modelling results
We are unaware of published work presenting lava flow modelling using a DSM. The closest equivalent is Charbonnier et al. (2018)'s study that models lahars. Charbonnier et al. (2018) found that when the lahar model FLO-2D is calibrated with geological data, a modelled lahar can be highly sensitive to sub-meter obstacles. We present lava flow modelling in an urban area and provide a quantitative measure of how the built environment alters the lava flow model results. There are locations where a DSM is unnecessary, and indeed its inclusion could needlessly increase processing times. For example, a DSM is unlikely to be required if the area modelled coincides with a park; indeed, a DSM of a forested park would include combustible trees which would be treated as permanent barriers, which likely does not reflect reality. In order to determine whether the results from a DSM or DEM are more appropriate to use given the wider hazard assessment context of the modelling, we propose using two measures to compare the DEM and DSM results. The first measure is a modified version of the Jaccard coefficient: where A DEM is the area covered by the lava flow when modelled on the DEM and A DSM is the area covered by the lava flow when modelled on the DSM. When the resulting overlap is high, then the DSM is not adding value to the simulations run on the DEM. Conversely, when the overlap is low, the implication is that human modifications to the land are resulting in a difference in the modelled lava flow areal footprint. Thus, the lower the percent overlap, the more important it is to consider what is causing the difference between the two footprints. If structures included in the DSM are likely to behave as barriers (e.g. concrete buildings), then the modeller should critically consider if modelling on a DSM is possible. The second measure compares the maximum run out lengths of the footprints: where L DSM is the runout length of the DSM results in the direction of the maximum runout length of the DEM results and L DEM is the maximum runout length of the DEM results. The lower the percentage is, the less overlap of the results there is in that direction.

Results
The lava flow model inputs are arranged by scenario in Table 2 while the model outputs are visually shown in Fig. 2 and compared in Table 3.

Discussion
Lava flow hazard models and computing power have greatly improved over recent decades, allowing for more detailed lava flow hazard assessments to be conducted on a variety of surface models.

Importance and implications of surface model selection
Although the selection of a lava flow method is important, choosing the correct elevation model for the project also greatly influences the accuracy of the model's outputs. When modelling hazards, it is standard to use a DEM (e.g. Barca et al. 1993;Scifoni et al. 2010;Favalli et al. 2012). Occasionally, the DEM may be modified to incorporate a manmade barrier to determine if the barrier could protect an area of significance (Fujita et al. 2008;Scifoni et al. 2010). Although using a modified DEM is not common practice, the topography in some volcanic areas has been extensively modified (e.g. Charbonnier et al. 2018;Hayes et al. 2018). Thus, DEMs may not be the most accurate representation of the obstacles a lava flow may encounter (Kereszturi et al. 2012), and a DSM may be more appropriate.
In our AVF case studies, we presented lava flow hazard footprints on both DEMs and DSMs. In each case, the hazard footprint generated on the DEM had a longer runout length (≥2% to 24% longer) than the corresponding DSM results (Fig. 2). The DSM results also frequently had a corresponding equal or thinner lava flow thicknesses the DEM results (Table 3). The shorter runout length and thinner footprints in the DSM results both have implications for evacuation zones and expected building damage. Shorter runout lengths may necessitate smaller evacuation or , so obstacles such as timber-framed buildings, represented in the DSM (Fig. 3), may be destroyed during the eruption and no longer affect the lava flow path (e.g. Williams and Moore 1983;Jenkins et al. 2017). Conversely, natural, topographic obstacles are unlikely to disappear during an eruption. Therefore, for lava flow modelling, a surface model that combines a DEM with aspects of a DSM (i.e. a DSM that does not include trees, or buildings where the primary building material has a low combustion point) may be the most accurate representation of a modified terrain. The original LAS dataset would be needed to develop such a hybrid DSM, and the endeavour is likely to be time consuming and to require familiarity with local construction practices. Ideally, lava flow models would allow for changes of the surface model from one iteration to the next to represent progressive damage to the built environment. This will require further research about lava flow impacts. Thus, we suggest that DSMs should be strongly considered when simulating lava flows advancing through heavily modified areas. In lightly modified areas and in areas where the dominant construction material has a low combustion point (i.e. timber-framed buildings), a DEM could still be more appropriate, especially if the final areal footprint is important.
Using a DSM can have important planning implications, but using a DSM also comes with a computational cost. DSMs are frequently produced using much higher resolution data than DEMS as surficial features tend to be on metre scales. Most lava flow models were not created with metre scale surface models in mind (Mossoux et al. 2016). Thus, processing times are often substantially longer when modelling on a high spatial resolution surface model than the low-resolution surface models traditionally used. Depending on the model used and computing power available, rapidly generating lava flow model outputs, as is required during a crisis, may not be possible on a DSM due to the higher resolution required of the topography model (or a DEM with equivalent spatial resolution). Despite these drawbacks, DSMs may be more appropriate input data in certain circumstances, such as when planning exercises.

Quantitative evaluation of surface model selection
Previous applications of the Jaccard coefficient in the context of hazard modelling has been used to assess how accurately a hazard model's output matches the scenario the model is attempting to replicate (Cordonnier et al. 2016;Mossoux et al. 2016;Favalli et al. 2012). We have adapted it to analyse the overlap between the DEM areal footprint and the DSM areal footprint. Following the lead of Cordonnier et al. (2016) and Dietterich et al. (2017), we do not provide guideline Jaccard coefficient percentages to determine when to use a DSM versus a DEM, as the appropriate cut-off percentage undoubtedly varies from one project to the next. The Jaccard coefficient indiscriminately calculates the overlap between two areas. If one is modelling on a DSM, it is important to note that the size of the kipuka (any area that is not inundated but is surrounded by a lava flow) in the results. If a building survives an eruption (forming an anthropogenic kipuka (e.g. the kipuka formed by a church during the eruption that began at Paricutin in 1943 (Luhr and Simkin 1993))) but is surrounded by lava, it will be difficult to access. Thus, some may consider such buildings to be useless. A kipuka can also be large enough to enclose multiple buildings, though (Neal et al. 2019). In such cases, the buildings may not sustain any damage although would still be difficult to access during and after the crisis. Both natural and man-made diversion structures could likely create such kipuka. Additionally, buildings surrounded by lava that survive an eruption and can still be entered may preserve smaller assets that can be evacuated when the lava flow is safe to cross. No matter the kipuka size, the kipuka will influence the Jaccard coefficient results (e.g. Fig. 2). For the aforementioned reasons, we feel that the Jaccard coefficient calculation should include the DSM areal footprint as is, i.e. with the kipuka created by the buildings. While including the diversion structure and building kipuka will decrease the total overlap compared to using only the perimeter of the areal footprints to calculate the Jaccard coefficient, the kipuka have planning implications. In some instances such as evacuation planning, the outer perimeter of the areal footprint is of more interest than the entire footprint.
Comparison of results with qualitative lava flow hazard footprints in Hayes et al. (2018) Most of the input parameters used by Hayes et al. (2018) are the same as those used here. In fact, Hayes et al. (2018) used a surface model, rather than an elevation model, to draw their qualitative lava flow footprints. The volume pulse is the only input parameter that is not shared between Hayes et al. (2018) and this paper. Despite very similar starting points, the resulting areal extents of the lava flows in Hayes et al. (2018) and this paper are quite different. First, the extents in Hayes et al. (2018) are notably longer (Fig. 4). This is likely because effusion rates, represented by the pulse volume parameter in MOLASSES, are not directly considered in Hayes et al. (2018). Discontinuous lava effusion was considered in Hayes et al. (2018) and represented by the stalling of the lava flow fronts; variations in effusion rates were, otherwise, not considered. In order to create MOLAS-SES footprints with similar extents to Hayes et al. (2018), a lower pulse volume would be necessary. Second, the surface models used in this paper had a higher spatial resolution than the surface model used by Hayes et al. (2018) enabling more precise topographic gradient calculations than were used in Hayes et al. (2018). Assuming accurate representation of the physical process, the paths lava followed in our case studies ( Fig. 2; Fig.  4b) are more realistic than those drawn in Hayes et al. (2018;Fig. 4a). Although both the effusion rate and spatial resolution discrepancies suggest that the areal extents presented here consider more nuanced data than those in Hayes et al. (2018), inputs to the MOLASSES modelling undertaken here could be improved. For example, we used the AVF lava flows' mean thickness, rather than the requested modal thickness, due to availability of data, but this assumes the AVF lava flow thickness distribution is Gaussian, which we deem unlikely. Thus, the lava flow modelling in this paper can be further updated as we learn more about AVF lava flows.
The case studies presented here provide more information (i.e. thickness maps on a metre scale) on the lava flows in the DEVORA scenarios than previously existed, so our model results can be used in a variety of preparedness activities including stakeholder workshops and evacuation and contingency modelling. The results of modelling conducted in this paper are aligned with Hayes et al. (2018)'s goal of using quantitative models for all of the considered hazards, although we acknowledge that qualitative methods may be more appropriate depending on the required outputs.

Conclusion
As settlements continue to encroach upon volcanoes, it is important to start or continue to prepare for future eruptions. Lava flow modelling has become increasingly sophisticated over the past few decades, but with the proliferation of models, selection of the most suitable model for a specific purpose becomes more difficult. Although most lava flow models have been developed to run on DEMs, DSMs capture all obstacles (irrespective of their physical properties) the flow may encounter. Thus, when lava flow hazard models are run on a DSM, the results will have considered the obstacles encountered by the flow. We illustrate the difference a DSM can make in three case studies in the AVF. In our case study scenarios, we found that lava flows modelled on a DEM generally have a longer runout length compared to lava flows modelled on a DSM. We use a modified Jaccard coefficient to quantify the difference between using a DEM and using a DSM. The most overlap (90%) occurred in the Ōtāhuhu scenario which is located in an area with minimal anthropogenic modifications. Thus, a DEM adequately portrays the topography in this scenario. Since there is much less overlap in terrestrial areas with generally high levels of urban development (e.g. the Birkenhead and Mt Eden scenarios), we would recommend using the DSM results in these cases. A similar decision would likely be made in most of the terrestrial Auckland Volcanic Field.