Thermal impacts of basaltic lava flows to buried infrastructure: workflow to determine the hazard

Lava flows can cause substantial physical damage to elements of the built environment. Often, lava flow impacts are assumed to be binary, i.e. cause complete damage if the lava flow and asset are in contact, or no damage if there is no direct contact. According to this paradigm, buried infrastructure would not be expected to sustain damage if a lava flow traverses the ground above. However, infrastructure managers (“stakeholders”) have expressed concern about potential lava flow damage to such assets. We present a workflow to assess the thermal hazard posed by lava flows to buried infrastructure. This workflow can be applied in a pre-defined scenario. The first step in this workflow is to select an appropriate lava flow model(s) and simulate the lava flow’s dimensions, or to measure an in situ lava flow’s dimensions. Next, stakeholders and the modellers collaborate to identify where the lava flow traverses buried network(s) of interest as well as the thermal operating conditions of these networks. Alternatively, instead of direct collaboration, this step could be done by overlaying the flow’s areal footprint on local infrastructure maps, and finding standard and maximum thermal operating conditions in the literature. After, the temperature of the lava flow at the intersection point(s) is modelled or extracted from the results of the first step. Fourth, the lava flow-substrate heat transfer is calculated. Finally, the heat transfer results are simplified based on the pre-identified thermal operating conditions. We illustrate how this workflow can be applied in an Auckland Volcanic Field (New Zealand) case study. Our case study demonstrates considerable heat is transferred from the hypothetical lava flow into the ground and that maximum operating temperatures for electric cables are exceeded within 1 week of the lava flow front’s arrival at the location of interest. An exceedance of maximum operating temperatures suggests that lava flows could cause thermal damage to buried infrastructure, although mitigation measures may be possible.


Introduction
Lava flow modelling has been undertaken for decades and, among other things, has been used to assess the hazard posed by flows to the built environment. For example, lava flow modelling during the 1991-3 Etna eruption (Dobran and Macedonio 1992) was used to assess whether the town of Zafferana Etnea could be inundated. Results indicated the town was at risk if a large enough volume of lava were erupted, prompting the design and construction of a lava flow diversion scheme (Barberi et al. 1993). Some lava flow modelling is conducted prior to an eruption to determine which areas could be threatened in a future or hypothetical eruption (e.g., Harris and Rowland 2001;Kereszturi et al. 2012;Mossoux et al. 2016;Gallant et al. 2018;Hayes et al. 2018;Cappello et al. 2019;Tsang et al. 2020). The results from modelling such lava flows allow for deeper analysis beyond areas likely to be inundated given the additional time available during times of quiescence. One use of lava flow modelling results prior to an eruption is in impact assessments (e.g. Deligne et al. 2015), which determine the level of damage particular exposed elements would experience in a future or hypothetical lava flow.
An impact is generally defined as the result of a hazard interacting with an exposed asset where the consequence is governed by the asset's vulnerability . Lava flows, especially those on shallow slopes, generally advance slowly enough that they do not pose a threat to life (Harris 2013;Kilburn 2015;Brown et al. 2017). However, immobile objects can be severely damaged by lava flows' high temperatures and pressures . In most impact assessments, lava flows are considered to cause complete damage if the flow comes in contact with an object (e.g. Wilson et al. 2014;Deligne et al. 2015;Hayes et al. 2018;Mossoux et al. 2019). However, anecdotes from the literature about global lava flow inundations, and from discussions with communities that have been impacted, suggest that destruction is not always complete (Tsang and Lindsay in press). For example, after the 1973 Eldfell eruption in Iceland, fish factories were able to be remediated and used again (Williams and Moore 1983). Similarly, a solid waste transfer station in Hawaii, USA was partially inundated in 2014 and returned to service in 2015 (Tsang et al. 2019a). Damage can also be caused when assets are not in direct contact with a lava flow. For example, lava flows have previously ignited forest fires that have damaged areas far from the flow itself (e.g. Suh et al. 2003;Davoine and Saint-Marc 2016;Staudacher et al. 2016), and caused up to a hundred fatalities when a petrol station exploded in 2002 (Brown et al. 2017). In yet another example, construction workers over 20 km downwind of the 2018 Lower East Rift Zone (Hawaii, USA) ocean entry commented that laze (i.e., the gaseous products when lava flows and ocean water interact) affected the safety procedures allowing them to work (e.g., whether they were allowed to work and what personal protective equipment was required) (Tsang and Lindsay 2019). Similarly, complex lava impact interactions are illustrated by fatalities during the 1930 Mt. Etna eruption when bystanders were killed by debris ejected by steam driven explosions when a lava flow overrode a submarine cistern (Francis and Oppenheimer 2004).
Stakeholders (e.g., infrastructure network and emergency managers) with assets near effusive volcanoes have expressed concern about the possibility of lava flows causing damage to buried infrastructure (Tsang and Lindsay in press;Tsang et al. 2019a, b) including electricity, fibre, potable water, stormwater, and wastewater (the 'three waters'), and gas networks. Stakeholders in Hawaii and New Zealand have identified that possible impacts to buried infrastructure are poorly understood despite the potential importance of such networks (Hawaiian Electric Light Company pers comm.; Hawai'i County Department of Water Supply pers. comm.; Auckland Council pers. comm.). This paper presents a workflow to assess the thermal hazard posed by lava to buried infrastructure. Damaged or operationally compromised infrastructure can result in service disruptions during and/or after an eruption. Only thermal hazards are considered here; mechanical hazards caused by other volcanic hazards, or by secondary lava flow hazards (including explosions caused by flow-water interactions and erosion by the lava flow) are not considered, although could produce impacts causing service disruption.

Lava flow hazard intensity metrics
A hazard intensity metric (HIM) is a characteristic of a hazard that can be measured and correlated to impact (Jenkins et al. 2014;Wilson et al. 2014). In the case of lava flows, the presence of lava (spatially defined by an areal footprint) is the most often used HIM as it is common to classify all assets in contact with the lava as fully damaged. Other HIMs commonly associated with lava flows are the pressure exerted by the flow and the temperatures to which the asset is exposed (e.g., Harris 2013). In the case of buried infrastructure, the gravitational basal pressure exerted by the lava flow would be well distributed over the large area covered by the basal crust, and thus, is unlikely to pose a substantial hazard. However, the heat transferred from the lava flow into the ground can be significant (e.g., Wilson 1962;Ishihara et al. 1989;Patrick et al. 2004;Fujita and Nagai 2016;Tsang et al. 2019b) and, therefore, could cause hazardous conditions for buried infrastructure. We note that thermal erosion of up to 10 cm/day has been documented (Peterson et al. 1994;Kauahikaua et al. 1998) and occurs when sufficient heat is transferred into the ground (e.g., Bussey et al. 1995Bussey et al. , 1997Fagents and Greeley 2001;Kerr 2009). However, as this process creates a mechanical, rather than thermal, hazard for buried infrastructure, it is out of scope for our study.

Thermal operating conditions of buried infrastructure networks
Ground heating is especially important for buried infrastructure as operating conditions are heat-sensitive. Buried infrastructure cannot operate at full capacity above certain temperatures (e.g., Rerak and Oclon 2017;Oclon et al. 2018;Eland Cables n.d.) nor at all above maximum operating temperatures (e.g., Evans 1981;Cinquemani et al. 1996;Schneider Electric 2002;Terpe 2017). Maximum standard operating temperatures can vary from ambient temperatures, in the case of three waters networks, to 100°C or above in the case of electric networks (e.g., Oclon et al. 2018;Eland Cables n.d.). While exceeding maximum standard operating temperatures is possible for some networks, there are disadvantages, such as decreasing the network's maximum possible load. The elevated substrate temperatures can also accelerate aging of buried infrastructure (e.g., de Leon et al. 2006;Bustamante et al. 2019). As buried infrastructure depths can vary from slightly above ground (e.g., when a gravity-driven pipe is exposed above ground; R. Roberts, pers. comm.) to tens of metres below ground (R. Joyce, pers. comm.), the amount of time between lava flow inundation and when critical temperatures are reached ranges substantially, from near instantaneous to weeks or longer (Tsang et al. 2019b). Stakeholders responsible for maintaining electricity networks located near volcanoes have expressed strong concerns about the amount and rate of heat transferred from the lava flow into the substrate as elevated temperatures can cause both temporary loss of function and/or long-term damage.
The workflow presented in this paper is designed to be jointly implemented by scientists and local stakeholders to minimise dissemination of sensitive asset data. The method creates thermal profiles that can be simplified to make them easy to interpret, and that facilitate and support science-based stakeholder decision-making. We demonstrate how this workflow can be applied in a deterministic case study in the monogenetic, basaltic Auckland Volcanic Field (AVF), which underlies the city of Auckland on the North Island of New Zealand (e.g., Hopkins et al. 2020 and references therein).

Method
Prior to implementing our method, a deterministic volcanic hazard scenario is required (e.g., Deligne et al. 2015Deligne et al. , 2017Hayes et al. 2018Hayes et al. , 2020. The scenario could be an ongoing eruption (in which case the necessary volcanological data relevant to the hazard may be estimated or measured; e.g., Ishihara et al. 1989, Dobran and Macedonio 1992, Barberi et al. 1993 or a deterministic scenario (in which hazard data can be prescribed; e.g., Deligne et al. 2015Deligne et al. , 2017Hayes et al. 2018).
Once a scenario has been defined or selected, lava flow hazard modelling or in situ measurements can be used to calculate or determine the lava flow's (potential) dimensions. Next, discussions with stakeholders determine which infrastructure may be exposed and provide data on the infrastructure's operating conditions. Then, the temperature at the intersection of the lava flow and the infrastructure network is calculated. In the penultimate step, the lava flow-substrate heat transfer is modelled. Finally, the results are optimised for usability. A pictorial overview of the method presented in this section is shown in Fig. 1; Table 1 summarises required data. The remainder of this section provides guidance and commentary while expanding on what the modeller-stakeholder team must do in each step.
Step 1: establish lava flow dimensions Our method requires the extent and thickness of the (potential) lava flow. This can be achieved by a variety of methods, including field measurements, expert opinion and/or quantitative modelling (e.g., Miyamoto and Sasaki 1997;Cappello et al. 2016;Fujita and Nagai 2016;Gallant et al. 2018;Harris et al. 2016;Kelfoun and Vallejo Vargas 2016;Hayes et al. 2018). If a quantitative model will be used, a lava flow model or ensemble of models must be selected. Numerous lava flow models exist, each with its own strengths and weaknesses. We direct you to  for advice on selecting an appropriate lava flow model and surface or elevation model ). If using the Tsang (2020) selection process, the required model outputs for this workflow are an areal footprint (for Step 2), lava flow thickness (for Step 2), lava flow temperature (for Step 3), advancement rate of the flow (for Step 4), and duration of the lava flow's activity (for Step 4). The advancement rate and duration of the flow may need to be approximated using analogue eruptions given that, as of mid-2020, there are no lava flow models that output these variables. Computational fluid dynamics models could potentially be used to simulate lava flows (Dietterich et al. 2017), although they are computationally expensive and would require substantial work to set up and verify on lava flows. Thus, we have not considered them here.
After a model or model ensemble is selected, the model providing the lava flow dimensions should be run. If a model ensemble composed of an areal footprint model and a thermo-rheological model will be used, the thermorheological model will be implemented in Step 3. Once the areal extent, thickness, advancement rate, and duration of the lava flow have been calculated, it is time to determine which infrastructure networks may be at risk.
Step 2: overlay infrastructure and determine operating conditions The lava flow hazard model results will identify where the lava flow will inundate, providing the spatial distribution HIM. One then needs to identify which infrastructure is exposed to or buried by lava in the given scenario. In this step, the areal footprint created in Step 1 is overlain with the local buried infrastructure network's location(s) to establish where infrastructure intersects the modelled flow path. The thickness of the areal footprint at the intersection point should then be noted for use in Steps 3 and 4. There may be multiple intersection points, all of which should be modelled in the following Steps 3 and 4. Once intersection points have been identified, the length of the path from the vent to the infrastructure's location can be measured. If the areal footprint is a flow field rather than a channel, an approximate length along the axis of the flow field may be used. Additionally, the slope of the identified path should be measured.
To determine how exposed buried infrastructure networks could be impacted by the modelled lava flow, knowledge of standard and maximum operating temperatures for each network is required. Therefore, Step 2 also determines temperature bounds to target the results to specific infrastructure. In some cases, both the standard and maximum operating temperatures could be the same. The best method to determine these temperatures is to discuss with the infrastructure provider, but the following bullet points help provide the basis of a network-specific discussion about potentially critical temperatures (i.e., standard and maximum operating temperatures) for buried infrastructure and some potential consequences if critical temperatures are surpassed.
Electricity: The upper limit of standard operating temperatures is commonly 100°C (R. Joyce, pers. comm.); electric networks can operate in warmer conditions, though (Eland Cables n.d.). Maximum operable temperatures are based on the materials the cables are made of, insulated with, and bedded in. Above these temperatures, outages and damage (such as melted casings) will likely occur. Fibre: The operating temperatures of fibre cables vary based on the materials the cable is made of. Coated optical cables and leaded glass fibre cables operate in environmental conditions of up to 200°C (J-Fiber n.d.) and 600°C (Fiberoptics Technology Inc. 2019), respectively; whereas plastic fibre cables operate at lower temperatures (Fiberoptics Technology Inc. 2019). Maximum operating temperatures can vary from 482°C for leaded glass fibre cables to 70°C for plastic fibre cables (Fiberoptics Technology Inc. 2019). Above these temperatures, outages and damage (such as melted casings) will likely occur. Potable water: Standard operating temperatures are close to ambient although operation can continue above ambient temperatures. Excessive heat could result in melted gaskets (Tsang et al. 2019a), Table 1 Table summarising necessary data Step Notes

Establish lava flow dimensions
The lava flow's dimensions can be measured if the lava flow is already emplaced. Otherwise, we recommend numerical modelling of the lava flow dimensions; required data vary depending on the model(s) selected. See Tsang (2020) for a model selection process. As of 2020, all models or model ensembles that produce the outputs necessary to implement this method require a topographic model, the vent location, and the effusion rate.

Overlay infrastructure and determine operating conditions
Required infrastructure data are the infrastructure locations (ideally as a GIS shapefile) and operating conditions. If infrastructure data are too sensitive to provide the scientist, the outputs of Step 1 can be provided to the stakeholder for them to indicate the location(s) of interest. Required data for operating conditions are: 1. The upper limit of the standard operating temperature of the infrastructure of interest, 2. The maximum operating temperature if the asset can operate above the standard operating temperature.

Model thermal profile of lava flow
The temperature of the lava flow at the location(s) of interest are necessary. Some numerical models output this information across the entire lava flow, in which case the local temperature should be extracted at this step. Otherwise, a lava flow thermal model (e.g., FLOWGO) should be implemented or a temperature estimate should be generated. See Tsang (2020) for more information.

Model lava flow-substrate heat transfer
Finite element analysis heat transfer modelling (Ansys APDL (https://www.ansys.com)) has previously been used to model lava flow-substrate heat transfer; see Tsang et al. (2019b) for more information.

Streamline results based on operating conditions
This step combines the outputs of steps 2 and 4 into a user-friendly format.
Tsang et al. Journal of Applied Volcanology (2020) 9:8 increased pressure due to boiling, and steam. These effects mostly occur at 100°C (maximum operable temperature). Storm water & wastewater: Both of these networks have similar considerations to the potable water networks. Gas: Standard operating temperatures of plastic gas pipelines in temperate regions is 49°C (Plastics Pipe Institute 2010). The maximum operating temperature of a plastic gas pipeline is 60°C (Plastics Pipe Institute 2010) above which the plastic will begin softening (Performance Pipe 2019).
The consequences described above are direct impacts caused by thermal hazard posed by lava flows. Secondary impacts such as explosions due to gas leaks may also occur. Although potentially severe, they are beyond the scope of this study.
Step 3: model thermal profile of lava flow After the intersection point(s) has been located, the flow's temperature at the intersection site should be determined. If the lava flow's core temperature was calculated in Step 1, then the core temperature at the intersection point should be extracted from the results obtained in Step 1. Alternatively, other approaches could be employed. One option would be to use the scenario parameters, the outputs from Step 1, and the topographical variables determined in Step 2 in a thermorheological model such as FLOWGO to determine the core temperature of the flow at the location of interest. Another method would involve using the heat budget of lava flows determined from analogue eruptions to estimate the flow's core temperature. Hon et al. (1994) describe one such relationship based on the first 10 years of the Pu'u 'Ō'ō eruption beginning in 1983, while Wooster et al. (1997) provide a different model based on the 1991-1993 Mt Etna eruption. The calculated temperature is a key input parameter for the heat transfer modelling conducted in Step 4.
Step 4: model lava flow-substrate heat transfer In this step, the heat transfer from the lava flow into the substrate below is modelled. Similar to the preceding steps, there are many ways to approach this step, from hand calculations to finite element analysis programs (e.g. , Jaeger 1959;Harris 2013;Rumpf et al. 2013;Turcotte and Schubert 2014;ANSYS 2017;Simpson 2017). For terrestrial applications, we recommend using the ANSYS APDL lava-substrate heat transfer model in Tsang et al. (2019b). The benefit of this ANSYS APDL model is that the conductance across the pāhoehoe lavasubstrate contact has been well defined and considers axisymmetric heat transfer in single-phase media. This means half a cross-section is "constructed" out of a single material with bulk material properties and then rotated around an axis of symmetry, in this case the y-axis, to create a three-dimensional model. As this is a singlephase model (meaning bulk material properties are used), phase changes, including the vaporisation of pore water, are not modelled. Convection and radiation are applied on the upper, external boundary of the model, while the bottom boundary is modelled as conductive into more dry soil. All other external boundaries are well-insulated; internally, heat transfer is modelled as conductive (Tsang et al. 2019b). APDL uses Galerkin's method to discretise the first law of thermodynamics to model convection and conduction (Madenci and Guven 2006;ANSYS 2017). We used APDL's radiosity solution method which employs the Newton-Raphson method to implement Stefan-Boltzmann's law to determine radiative fluxes (ANSYS 2017). Using the axisymmetric heat transfer model described in Tsang et al. (2019b) requires quantifying the following: the thickness of the flow, the core temperature of the flow, the flow front's advancement rate, the duration of the activity, and the substrate characteristics (including moisture content). The flow's thickness and the temperature at the intersection point are determined in Steps 1 and 3, respectively. The substrate can be identified using a variety of measures including geological maps, field work, geotechnical reports, and satellite imagery. It is also important to consider the water table depth to estimate the moisture content of the substrate. The flow front advancement rate is used to determine how long it would take for the lava flow to reach the location of interest. Finally, how long lava is supplied to the intersection location needs to be determined (duration of "actively flowing phase" in Tsang et al. (2019b)). Few lava flow models calculate this, so a first-order estimate can be generated using the total eruption duration minus the amount of time it takes the flow front to reach the location. The "cooling-only phase" (Tsang et al. 2019b) starts when active lava is no longer supplied and ends when the soil is cooling and/or all of the lava has cooled. Determining the cooling-only duration may take several iterations of modelling. Should an extra-terrestrial application arise, the heat transfer model in Rumpf et al. (2013) may be most appropriate. Other finite element analysis options include other commercial packages such as COMSOL and Abaquas, the FEATool Multiphysics toolbox (https://www.featool.com/), or the method outlined in Chapter 8 of Simpson (2017). Alternatively, Harris (2013) describes how to calculate by hand the heat budget of a cooling lava flow including the temperature of the flow's basal crust. This temperature could then be applied in the contact temperature equation described in Jaeger (1959). Conduction in the substrate can then be computed.
While previous calculations of lava flows transferring heat into the ground could potentially be applied here, it is unlikely the specific situation being considered has already been modelled. Substrate properties vary substantially from one location to another, as do road and footpath construction. Thus, we strongly recommend modelling the heat transfer in the specific location of interest rather than relying on previous results.
The temperature profile under the lava flow generated in this step provides more data than decision-makers need. The standard and maximum operating temperatures determined in Step 2 will aid in streamlining the results.
Step 5: streamline results based on operating conditions Once standard and maximum operable temperatures are known, the results from Step 4 can be adapted to a usable format for decision making. Given that temperature profiles change during the eruption, we recommend providing stakeholders a temporal series of profiles identifying four thermal zones of interest (see Fig. 4c for an example diagram). The first zone is likely deepest (i.e. farthest from the lava) and represents temperatures at or below the standard operating temperature (i.e., the coolest section of substrate, and outside the thermal hazard zone). The second zone (the moderate hazard zone) represents temperatures between the standard operating temperature and the maximum operating temperature. The third zone (the high hazard zone) represents temperatures above the maximum operating temperature and below the substrate's bulk solidus (if known). The fourth and final zone is the top substrate zone and represents temperatures above the substrate's solidus. This zone corresponds to the thermal conditions suitable for thermal erosion, although thermal erosion is also reliant on other factors and, therefore, may not occur (as is the case in Fig. 4c). It is key to include a legend summarising the expected conditions in each zone.
Lava flow inundations and potential mitigation measures will vary from volcano to volcano and from one event to the next. This versatile workflow can be adapted to determine the heat transfer under a basaltic lava flow in a wide variety of substrates.

Results: Auckland Volcanic Field case study
We demonstrate how this workflow can be applied in a deterministic case study in the monogenetic, basaltic AVF (Fig. 2), which underlies the city of Auckland in the North Island of New Zealand. The city, built since the last AVF eruption, is home to 1.6 million people (Stats NZ 2019a) and generates a third of the country's gross domestic product (Stats NZ 2019b). Thus, both Auckland and New Zealand could be severely impacted by a future eruption (e.g., Deligne et al. 2015Deligne et al. , 2017Blake et al. 2017;Hayes et al. 2017;McDonald et al. 2017). The Determining Volcanic Risk in Auckland (DEVORA) research programme (http://www.devora.org.nz) was launched in 2008 to study the volcanic risk posed to the city by proximal and distal volcanoes. One research output has been eight hypothetical, multi-hazard, deterministic eruption scenarios (Hayes et al. 2018), developed following stakeholder demand ).

DEVORA scenarios
The DEVORA scenarios are eruptive sequences that represent the full range of possible eruptive phenomena and hazard intensities that could be expected during an AVF eruption (Hayes et al. 2018. They were developed by the DEVORA research team through multi-hazard modelling and a series of workshops. Each scenario is composed of multiple hazards based on the hazard's frequency in past AVF eruptions and environmental conditions (e.g., groundwater depth; Hayes et al. 2018). We selected the Birkenhead scenario for our case study illustrating the initial stages (i.e., first lava flow of the scenario) of disruption, as the lava flows in this scenario could potentially impact major buried infrastructure (i.e., national electrical transmission lines). For more details, see Hayes et al. (2020). For variables not defined in Hayes et al. (2018), data sources are indicated in Table 2.
Step 1: establish lava flow dimensions We followed the outputs-focused framework for selecting lava flow models presented in  and identified that MOLASSES (Gallant et al. 2018) and FLOWGO (e.g., Harris et al. 2016;Chevrel et al. 2018) are best suited to our purposes. We used MOLASSES, a digital surface model from Land and Information New Zealand, and the values in Table 2 to model the areal footprint and thickness of the first lava flow in the Birkenhead scenario.
The resulting lava flow footprint is shown in Fig. 3 and shows the lava flow advancing to the north and east of the vent. The easternmost extent of the lava flow advances into the inter-tidal zone of the Waitemata Harbour (Fig. 2). This could result in the formation of a lava delta (e.g., Skilling 2016) and/or explosions caused by lava-water interactions (e.g., Mattox and Mangan 1997;Hamilton et al. 2010;Noguchi et al. 2016). Lava deltas and lava-water interactions could pose mechanical hazards to surrounding locations although are not discussed further in this paper.
Step 2: overlay infrastructure and determine operating conditions After modelling the expected areal footprint, we used open-access data (Transpower New Zealand 2019; Watercare Services Ltd. 2019) to overlay the potable water pipes, wastewater pipes, and electricity transmission cables in the area. Local government data indicated that there are no major (i.e., national transmission) gas pipelines or fibre cables in the area. Only the electric transmission cable, which is located 250 m from the vent, is part of a national grid, so for illustrative purposes the remainder of this paper will focus on this transmission electricity cable. There is no local redundancy for the electricity transmission cable, which runs under the highway on the right side of Fig. 3. The average slope between the vent and the electric cable location varies between 1.4 and 1.9°based on the final elevations calculated by MOLASSES (i.e., the slope of the lava flow as simulated by MOLASSES). Above the electric cable, the modelled lava flow is 17.5-m-thick and 54-m-wide. Given that MOLASSES produces a final areal footprint without any temporal dimension, we rely on the advancement rates and active durations published in Hayes et al. (2018), which were based on expert judgement. The upper range of the standard operating temperature of the electric cable is 100°C (Transpower New Zealand pers. comm.). Some electric cables can continue operating at a lower capacity above 100°C. For illustrative purposes, we will use 150°C as the maximum operating temperature, which is the maximum operating temperature of silica cables (Eland Cables n.d.).
Step 3: model thermal profile of lava flow Because MOLASSES does not calculate lava flow temperatures, in this step we model the temperature of the lava flow using FLOWGO and the values in Table 3 along the path described in Step 2.
FLOWGO calculates the temperature of the core of the lava flow at our location of interest to be 1199°C.
Step 4: model lava flow-substrate heat transfer Using the outputs of the previous steps and the values in Table 4, we run the heat transfer model described in Tsang et al. (2019b) to simulate the heat transfer from the lava flow into dry soil. We used ANSYS APDL for the heat transfer modelling at local stakeholders' request (e.g., Transpower engineers). Due to the proximity of the location of interest to the harbour, the water table is likely close to the surface, but data from boreholes in the area indicate that the water table is more than 5 m deep at the location of interest (New Zealand Geotechnical Database 2020). Thus, we have assumed that the water table depth is more than 5 m deep, and modelled the bottom boundary of our substrate as if it were resting on top of dry soil, rather than on water-saturated soil. This decision is further supported by the proximity of the lava flows to the vent, suggesting the ground may have been desiccated by the eruption, similar to previous AVF eruptions ( Kereszturi et al. 2014b). Thus, we have opted to model the lava flow on top of dry soil. Since the "soil" used in Tsang et al. (2019b) is based on geotechnical reports of soils in Auckland (Taihan New Zealand Ltd. 2010; Rifareal 2011), we use the soil parameters specified in that publication without alterations. Variations in soil conditions with increasing depth are not considered in this example but could be if such data were available. Hayes et al. (2018) provide a flow front advancement rate of 30 m/hr., and state that lava effuses from the vent for 3 days. Assuming the flow does not stall, the flow front should reach our location of interest in 8 h and 20 min, which we rounded to 8 h to make our flowing phase duration estimate conservative. The flowing phase at our location of interest is, therefore, 64 h. The cooling-only phase duration was eventually set at 9 months based on when heat transfer modelling indicated that substantial heat was no longer being transferred from the lava to the substrate. (See Tsang et al. (2019b) for more information about first order controls on substrate thermal profiles).
A temporal series of temperature profiles is shown in Fig. 4 (full results are shown in Online Resource A). By  the end of the emplacement of the flow and the end of the active flowing phase, the lava flow has only heated the upper metre of substrate; the entire substrate column is still under 100°C. A week after the lava flow emplacement ends, the top 2.5 m of substrate has been heated above ambient temperatures while the upper 1.5 m is above 150°C. As the lava flow cools, it continues to heat the substrate. One month after the emplacement ends, the lava has heated the top 3.5 m of substrate above 150°C. Given that buried infrastructure is commonly shallow (less than 3.5 m deep), this indicates most buried infrastructure will be surrounded by substrates above maximum operating temperatures within a month of lava's emplacement. By the end of 6 months, the entire 5 m substrate column has been heated to above 150°. The heating continues past 6 months. Nine months after the active flowing phase ended, the entire substrate column is still above maximum operating conditions with the upper few centimetres heated above 800°C.
We do not provide a temperature for longer durations because the lava's maximum temperature at 9 months (found in the basal crustal area) is 200°C lower than its original emplacement temperature and the temperature in the upper 5 cm of soil has changed less than a degree in the final three simulated months. The simulated high temperatures are comparable to those found in thick lava flows after an eruption; for example, the Eldfell lava flow field was still cooling 10 years after the eruption ended (Williams and Moore 1983). The small temperature change in the upper soil during the last 3 months indicates that heat is being transferred into the soil at a slower rate than in the first few days after the eruption ended. Although soil temperatures could continue to increase after 9 months, the modelled temperatures in the top 5 m of substrate are already more than four times maximum operating conditions, meaning the buried infrastructure networks have not functioned in months. The mesh size (Table 4) is an important parameter in finite element analyses. While an increased mesh size allows for more rapid calculations, it also increases the resulting error. Given that we have increased the mesh size compared to the model presented in Tsang et al. (2019b), we have calculated a grid convergence index (GCI), extrapolated error, and approximated error. These uncertainties represent the error introduced by turning the partial differential equations into discrete systems of equations and are 1.1 × 10 − 14 , 8.5 × 10 − 15 , and 0.019%, respectively.
Step 5: streamline results based on operating conditions In the last step, we adapt the temperature profiles we created in Step 4 with the thermal bounds identified in Step 2. We have taken the first four profiles in Fig. 4b to create Fig. 4c, in which standard operating temperatures are shown in grey, operable temperatures above standard are shown in orange, and temperatures too hot to continue operations are in red. In this case study, the substrate temperatures are not high enough to initiate thermal erosion, so this zone is not displayed or included in the legend. The resulting visual depictions of the thermal hazard under the lava flow can aid in decision making.

Discussion
We have outlined a method to evaluate the thermal hazard posed to buried infrastructure by lava flows in a deterministic scenario using lava flow hazard models, heat transfer modelling, and stakeholder input. We then applied the model to the first lava flow in the Birkenhead scenario developed for the Auckland Volcanic Field by Hayes et al. (2018). We found that within a month, most buried infrastructure in the Birkenhead residential area will be exposed to well above operable temperatures. The substrate surrounding the electricity transmission cable, assumed to be buried within the top 1.5 m of substrate (Transpower pers. comm.) will be above maximum operating temperatures in less than a month (Fig. 4b).

Considerations when implementing this method
This method is only as accurate as the data and models used in it. While some data can be easily measured or are commonly collected, not all of the input variables (e.g., variables listed in Tables 2, 3, and 4) may be known. In many cases, both the lava flow and the asset characteristic datasets are incomplete. Thus, data for generic basaltic flows, from analogue volcanoes and eruptions (e.g., emissivity in Table 3) and/or expert judgement may be required. This introduces a source of uncertainty. Depending on when the generic or analogue data is used will determine how the uncertainty will propagate. Additionally, the uncertainty will propagate differently depending on which models are used.

Our model selection
In our case study, we used MOLASSES to model the areal extent, FLOWGO to model lava temperature, and ANSYS APDL to model the heat transfer, following Tsang et al. (2019b) and Tsang et al. (2020). We briefly comment on model uncertainty here. We selected MOLASSES and FLOWGO according to the process described in Tsang. (2020); both models have been Table 4 Input values for the lava-substrate heat transfer modelling. The mesh size was doubled compared to the mesh size in Tsang et al. (2019b) given the thicker lava flow here. The lava temperature is based on Step 3, and the ambient temperature is from CliFlo (a national New Zealand climate database; https://cliflo.niwa.co.nz/). Modified from Tsang et al. (2019b)  These parameters were not treated as temperature dependent or only a bulk coefficient was used. Thus, a single value at all temperatures was used d This is the core temperature of the lava flow Fig. 4 a Schematic diagram of the process being modelled in Step 4. b Mean temperature profiles in the 5 m of dry soil below the Birkenhead case study lava flow at 4 days (i.e., the end of the actively flowing phase), and at 1 week, 1 month, 6 months and 9 months since the start of the cooling-only phase. c Results shown in (b) presented in the context of the operating temperatures of interest at the end of the actively flowing phase (4 days), and a week, a month, and 6 months after the start of the cooling-only phase. Nine months after the start of the cooling-only phase is not shown as it is the same as for 6 months (i.e., also a solid red bar). A schematic diagram along the lines of (c) can be provided to stakeholders to aid in decision making tested in a variety of conditions. MOLASSES has been validated on the 2012-2013 Tolbachik (Kamchatka, Russia) eruption (Kubanek et al. 2015). Additionally, Dietterich et al. (2017) benchmarked MOLASSES against four other computational fluid dynamics codes and found that MOLASSES was one of the top two performing codes in four of the eight benchmark tests. Less benchmarking has been undertaken on FLOWGO, but it has been verified on five volcanoes around the world in 12 eruptions (Harris et al. 2007Rowland et al. 2005;Harris andRowland 2001, 2015). A sensitivity analysis conducted on two Mt Cameroon eruptions is available in Wantim et al. (2013). ANSYS APDL is a commercial finite element analysis software package that has been extensively verified and validated. (For more details about ANSYS APDL, see Online Resource 2.) Tsang et al. (2019b) presents an additional verification case for our heat transfer scenario.
In their implementation, they found a GCI of 0.015%, an approximate relative error of 0.220%, and an extrapolated error of 0.116%. Our implementation is slightly different due to an increased mesh size, but our GCI, approximate relative and extrapolated errors are even smaller still. Thus, ANSYS APDL is not introducing noteworthy error. Our heat transfer results indicate substantial heating to the substrate below; these results are in-line with previous work. For example, Bogue and Glen (2010) studied two stacked lava flows and found evidence that lava 3 m below the contact between the flows was heated to temperatures in excess of 300°C as the upper flow cooled. This was supported by one-dimensional conduction calculations. They concluded that the lower flow must not have been fully cool for this situation to occur. Similarly, Manley (1992) and Forman et al. (1994) model thermal diffusion and conduction, respectively, and find that substantial heating of underlying sediments can occur. While none of these findings are directly comparable due to different assumptions and notable differences in substrate properties, they support that substrates beneath lava flows can be heated considerably above ambient temperatures.

Lava flow morphotypes
An important note is that the case study illustrating this method used lava flow models written for multiple morphotypes. Gallant et al. (2018) do not specify what type of lava flow MOLASSES models although it has been applied to the Eastern Snake River Plain volcanic province (USA), which is characterised by large volume flows advancing across gentle slopes (Gallant et al. 2018) and to the 2012-2013 Tolbachik eruption which was a small volume flow on steeper slopes (Kubanek et al. 2015). In the benchmark tests conducted by Dietterich et al. (2017) on computational fluid dynamics models including MOLASSES, analogue molten rock flows which simulate pāhoehoe morphologies were also modelled with good agreement. MOLASSES can therefore potentially be applied to a variety of situations and types of flows. FLOWGO, conversely, can only be applied to channelised a 'ā' flows (e.g., Harris et al. 2016;Chevrel et al. 2018). The ANSYS APDL modelling presented in Tsang et al. (2019b) has only been verified on analogue pāhoehoe experiments and under a pāhoehoe flow. Although FLOWGO and ANSYS APDL have each been verified on different flow types, measurements from under lava flows indicate that the lava flow morphotype does not profoundly influence the decrease in lava flow temperature close to the vent (Glaze and Baloga 2016). Hon et al. (1993Hon et al. ( , 1994 measured the core temperature of pāhoehoe Kīlauea flows and determined that the flow's core temperature tends to drop 1°C per 10 km (Hon et al. 1994). At a distance of 250 m from the vent, this translates to a 0.025°C decrease in the lava flow's core temperature. FLOWGO (which is used to model lava flow channels) calculates a temperature drop of 1°C 250 m from the vent. This demonstrates that, regardless of which method is used to determine the temperature drop, the lava flow's core temperature is negligibly cooled in the first quarter kilometre from the vent. Therefore, given the core temperature estimates are within one degree of each other regardless of the lava flow type, FLOWGO results can be used with the heat transfer model. This is the first implementation of a workflow to consider the lava flow hazard to buried infrastructure using verified models. Although there are limitations as described above, the caveats do not seem to greatly affect the final result, i.e. the zones produced by Step 5.

Mitigation measures for buried electricity networks
There are several potential mitigation options that may decrease the vulnerability of buried electricity infrastructure to lava flow thermal hazard. We acknowledge that other volcanic hazards (including precursory earthquakes and rootless vents in lava flows) in the eruption sequence could also damage buried infrastructure; potential impacts posed by other volcanic hazards and their associated hazard intensity metrics should be critically considered prior to implementing a mitigation option. The following mitigation options are based on discussions with stakeholders and published literature, including a review of previous lava flow inundations (Tsang and Lindsay in press).
Some mitigation measures can be applied prior to an eruption to prolong network functionality. Discussions with electricity providers in Hawaii and New Zealand have prompted suggestions such as: Burying cables deeper (Klimenta et al. 2018): Burying cables deeper is most easily done when a cable is initially laid, but it increases the installation cost and makes maintenance more difficult (Salata et al. 2015). Wrapping cables in thicker or different insulation (Terpe 2017; Eland Cables n.d.): Alternative insulation and bedding materials can be placed during installation. A downside is that making the cables more insulated will negatively affect the cables' standard operations (de Leon et al. 2006;Salata et al. 2015;Bustamante et al. 2019). Changing the bedding of the cables (Ocłoń et al. 2015;2016): Similar to insulation, the bedding of the cables (i.e., the materials in which the cables are buried) can influence how the cables heat and cool. Using different types of cables: Other cables with different temperature ratings could be used (Cinquemani et al. 1996).
Considering the siting of buried networks: The siting of the cables is important as cables in wet soil will likely heat above operable conditions faster due to convection cells (Baker et al. 2015;Tsang et al. 2019b), suggesting it would be preferable to site cables in areas where the groundwater table is deep. Additionally, siting networks under footpaths and roads could be a potential mitigation measure as the analogue experiments in Tsang et al. (2019b) suggest the heat dissipates across associated bedding layers. The temperature profiles under molten rock on a road profile measured by Tsang et al. (2019b) suggest that substantial heat from the lava flow is initially used to cause a phase transition in the asphalt of the road, rather than simply heat the substrate. Eventually, the temperatures followed a similar trend to those without the road, but burying cables under roads may prolong how long the cables can function, assuming they are situated sufficiently deep below the asphalt layer. The road bedding will also impede the heating of the soil surrounding the electric cable. This has the added benefit of being relatively easy to restore after the eruption, given that construction contractors in Hawaii have stated that roads can be easy to remediate, especially if they have been sealed, as the lava and sealant layer easily peel off the road bedding (Tsang and Lindsay 2019). Build redundant systems: Network redundancies allow for damage to one section of the network without compromising the entire network.
In summary, initial installation conditions can have considerable influence on how resilient the network may be to lava flow impacts. While these considerations focus on actions that could be conducted by electricity companies, most of them may also be applicable to other types of buried infrastructure networks.
Assuming mitigation measures have not already been implemented prior to an eruption, other mitigation measures can also be considered. Some of the following suggestions would aim to prolong the electricity cables' operations while others would involve rapidly building redundant systems: The cables' load could be decreased (Johnson and Propst 1989) and/or the network could be cycled on and off, causing brownouts (i.e., rolling blackouts): Since many systems are run at the upper limit of their standard operating temperatures to maximise their loads, such actions would decrease the heat being released by the electric cable or temporarily remove a heat source, respectively. Although this is not ideal, it would enable the network to continue functioning for longer. By minimising the heating of the asset, sheathing and/or gaskets may take longer to melt after lava traverses the ground above, increasing the potential the asset may survive the eruption even if it cannot function the entire eruption. This could be especially effective in lava flow events in which the lava flow is short-lived and/or thin, since less heat will be transferred into the ground in such cases. Once the substrate temperature has been elevated above maximum operating temperatures, little can be done to maintain service and it is questionable whether the cable can be returned to service. Build a redundant electric transmission line (i.e., aerial electrical network): A electricity line maintained full functionality while the ground below it was inundated by lava in Hawaii in 2014 (Tsang et al. 2019a). The continuity of operations has been attributed to the power pole protection measures that were built, the flow's pāhoehoe morphology, and its final width. The power pole protection measure increased how long the power pole remained standing once the lava flow arrived at the location of interest. The flow's cohesive pāhoehoe upper crust meant less heat was released into the air above the flow (i.e., the air surrounding the electric line) compared to an 'a'ā flow. Finally, the flow's width meant fewer power poles were required to suspend the line above the lava. Thus, if an area is at high risk of being inundated by lava and a transmission line or cable must be located in that area, overhead electric lines may function longer than buried electric cables will, assuming the flows are narrow and pāhoehoe in nature, and that subaerial networks are not impacted by other volcanic hazards.
Many of these suggested mitigative measures require forethought to ensure proper materials are at hand, so it is important to consider the limitations of the method presented prior to implementation.
In the Birkenhead scenario, lava is actively supplied to the location of interest for less than three days before the start of the cooling-only phase. Assuming that that the transmission cable is buried at 1.5 m, a standard burial depth in Auckland (Transpower pers. comm.), three days after the cooling-only phase commences (a week after the flow front reaches the cable's location), the temperature at 1.5 m is already over 100°C (Online Resource A). Since this cable is already installed, cycling the network may help maintain service through the first week of the hypothetical eruption. Within the first month of the lava flow cooling, the transmission cable will be surrounded by or in substrates above operable temperatures. At this point, service will have to be provided via an alternate route.

Future work
While our method enables the evaluation of the lava flow thermal hazard to buried infrastructure, more work can be done to improve this method and make it more broadly applicable. For example, this method relies on a heat transfer model with several assumptions (i.e., a pāhoehoe morphology, static thickness, no horizontal heat transfer, and no phase changes) that could be refined (Tsang et al. 2019b) and on lava flow footprint models that simulate different types of lava flows. Assumptions in the heat transfer modelling could be refined with additional field work to understand 'a'ā lava flows, with more laboratory lava flows, and/or with further heat transfer modelling. The lava flow footprint modelling could be improved if models to simulate pāhoehoe lava flows were created. All of these models could be further developed to improve hazard assessments. Additionally, the method introduced here could be automated and packaged in a graphical user interface. Doing so would allow stakeholders to undertake their own hazard analysis, which can be used in business continuity and emergency management planning. As lava flow models are developed and modified to output additional hazard intensity metrics, such as pressure exerted by the lava flow on the substrate, the method could also be edited to provide a more holistic hazard assessment for assets in the built environment. Buried infrastructure impacts caused by other volcanic hazards could also be considered with additional research.

Conclusion
In summary, we have presented a method to assess the thermal hazard posed by lava flows to buried infrastructure. This method draws on knowledge about the local infrastructure substrates and networks, areal and thermorheological lava flow models, and heat transfer modelling.
We displayed how the method can be applied in a case study in Auckland, New Zealand using the Birkenhead DEVORA scenario. The lava flow areal footprint created by MOLASSES revealed that the national electricity transmission cable in the area could be affected in this scenario, and subsequent heat transfer modelling determined that if the cable is buried 1.5-m deep, the cable would only continue operating between a week and a month after the eruption onset. Mitigation measures could be employed to prolong possible operations although have not yet been tested. This method can be employed not only by stakeholders to inform decision-making but can also be used to simulate potential mitigation measures.

Additional file 2.
Abbreviations AVF: Auckland Volcanic Field; DEVORA: Determining Volcanic Risk in Auckland; GCI: Grid convergence index; HIM: Hazard intensity metric