Probabilistic hazard modelling of rain-triggered lahars

Abstract

Probabilistic quantification of lahar hazard is an important component of lahar risk assessment and mitigation. Here we propose a new approach to probabilistic lahar hazard assessment through coupling a lahar susceptibility model with a shallow-layer lahar flow model. Initial lahar volumes and their probabilities are quantified using the lahar susceptibility model which establishes a relationship between the volume of mobilised sediment and exceedance probabilities from rainfall intensity-frequency-duration curves. Rainfall-triggered lahar hazard zones can then be delineated probabilistically by using the mobilised volumes as an input into lahar flow models. While the applicability of this model is limited to rain-triggered lahars, this approach is able to reduce the reliance on historic and empirical estimates of lahar hazard and creates an opportunity for the generation of purely quantitative probabilistic lahar hazard maps. The new approach is demonstrated through the generation of probabilistic hazard maps for lahars originating from the Mangatoetoenui Glacier, Ruapehu volcano, New Zealand.

Lahars are among the most hazardous volcanic processes and are responsible for a large proportion of volcanic fatalities (Auker et al. 2013). Reliable information on the likelihood of lahar occurrence and the resulting inundation area is critical for the mitigation of risks posed by lahars (Pierson et al. 2014). Lahar hazard assessments typically provide this information in the form of hazard maps quantifying the probability and extent of potential lahars to varying degrees (Calder et al. 2015). In these assessments, models that approximate lahar behaviour and/or run-out are frequently used (e.g. Aguilera et al. 2004; Carranza and Castro 2006; Darnell et al. 2013; Thouret et al. 2013; Pistolesi et al. 2014; Córdoba et al. 2015). Commonly used models include Laharz (Iverson et al. 1998; Schilling 2014), an empirical model relating lahar volume to cross-sectional and planimetric inundation area; the two-phase shallow layer model of Pitman and Le (2005) implemented in the Titan2D toolkit (Pitman et al. 2003; Patra et al. 2005); and single phase rheology approaches such as Flo-2D and Delft3D, used in Caballero et al. (2006) and Carrivick et al. (2009, 2010). Regardless of methodology, the accuracy and output of these methods is dependent on the value and accuracy of model inputs. The initial lahar volume or volumetric flux is one such input that determines initial lahar size and gravitational potential energy. Since lahar size and energy transfer are important factors controlling lahar behaviour and runout (Lube et al. 2012), an accurate estimate of initial volume is therefore crucial for accurate and reliable hazard footprint estimates. Commonly, initial lahar volumes are estimated from volumes of previous lahar events or rely on expert judgement. However, historic data is often incomplete or can be irrelevant under different environmental conditions which limits probabilistic lahar hazard estimates.

Recently, Mead et al. (2016) presented a physically based model for determining initial volumes of rain-triggered lahars. Lahar susceptibility, defined as the probability of an initial lahar volume at a specific location, was determined through assigning annual probabilities to lahar volumes using rainfall intensity-frequency-duration (IFD) curves. The potential of coupling lahar susceptibility model outputs to lahar flow model inputs in order to quantify lahar hazard probabilistically was discussed, but not demonstrated. Here, the work of Mead et al. (2016) is extended by coupling lahar susceptibility outputs to the two-phase shallow layer model of Pitman and Le (2005). Outputs of the combined model in the form of probabilistic hazard maps are demonstrated for lahar scenarios originating from Mangatoetoenui Glacier, Mt. Ruapehu, New Zealand. This work is intended to demonstrate an approach for probabilistic estimates of lahar hazard and highlight areas of research needed to enhance the reliability of these hazard estimates.

Methods

The probabilistic lahar hazard approach proposed here relies on the combination of two modelling methods. Lahars are simulated as a two-phase flow using Titan2D with initial volumes determined using the lahar susceptibility model of Mead et al. (2016). Details of each method are summarised below; for full details readers are referred to Patra et al. (2005) and Pitman and Le (2005) for the Titan2D toolkit and Mead et al. (2016) for the rain-triggered lahar susceptibility model.

Rain-triggered lahar susceptibility

In the rain-triggered lahar susceptibility approach described in Mead et al. (2016), lahar volumes are estimated from deposit properties, rainfall intensity and rainfall duration using a combined shallow landslide and overland erosion model. The total initial volume of the lahar is calculated in the combined model as the sum of the material mobilised through overland erosion and shallow landsliding for a given rainfall intensity and duration. In the overland erosion model, the height and motion of infiltration excess rainfall (i.e. net rainfall after infiltration into the deposit) is simulated using a depth averaged shallow water (SW) approximation. The entrainment of volcanic sediment is calculated from the height and velocity of the overland flow using the Meyer-Peter and Müller (MPM) bedload transport model (see Castro Díaz et al. (2009)), and the total volume of mobilised volcanic sediment is calculated at the end of the rainfall duration. The volume of volcanic material mobilised through shallow landsliding is calculated using the approach developed by Iverson (2000) where shallow slope failures are assumed to occur when gravitational forces are greater or equal to the resisting Coulomb stresses. In the Iverson (2000) model, resisting stresses are reduced in proportion to the rainfall infiltration rate (i.e. deposit permeability) and hydraulic diffusivity, which controls the transmission of pore pressure through the deposit. The depth, and consequently volume, of the shallow failure is calculated at the end of the rainfall duration. The probability of specific lahar volumes can be estimated using this method through the use of rainfall intensity-frequency-duration (IFD) tables. These tables express the probability of rainfall intensities occurring over a given duration and can be used to estimate the occurrence probability of specific lahar volumes being triggered by rainfall. These lahar volumes, when used as an input in Titan2D, can then be used to express lahar inundation in a probabilistic manner. Here, the mobilised volume estimates and recurrence intervals calculated in Mead et al. (2016) are used as the probabilistic input for the lahar flow model.

Lahar flow modelling

The choice of lahar flow model can have a large effect on the quality and reliability of probabilistic hazard estimates and therefore needs to be considered in conjunction with the needs of the hazard assessment. The prediction of lahar flow is difficult due to complex physical processes such as entrainment and deposition of sediment, changing solid concentrations and flow transformations (e.g. Doyle et al. 2009). While empirical models such as Laharz can be used to estimate lahar inundation areas (Schilling 2014), numerical models are able to provide more information (e.g. velocities and pressures) useful for estimating the intensity of the hazard (McDougall 2016). A range of numerical models have been used or proposed (e.g. Pudasaini 2012) for lahar flow modelling; however, the selection of numerical models and their parameters is complicated by the lack of universal constitutive laws governing lahar flow behaviour (McDougall 2016).

While input parameters for simple models that rely on bulk measures of the lahar (e.g. internal and basal friction, solid concentration) can be calibrated to previous lahar flows, in practice, the complexity of lahar events mean the calibration can only be optimised on one characteristic of the lahar to the detriment of others. For example, a reduction of friction values to match long runout lahars will result in excessively mobile lahars with high momentum in proximal (near source) zones. Recent, more advanced models of two-phase flow (e.g. Pudasaini 2012; Iverson and George 2014; Iverson and George 2016; Mergili et al. 2017) may improve predictions of lahar flow. However, these models rely on a large array of input parameters that are often not known or measurable a priori and, as yet, do not demonstrate the ability to model transformations in flow behaviour.

Probabilistic lahar hazard estimates are obtained through the use of volume and location outputs from the lahar susceptibility method as the input for lahar flow models. Consequently, the flow modelling approach can be chosen according to an assessment of modelling capabilities needed to accurately simulate the lahar hazard. To demonstrate this approach, we chose to use the two-phase material model of Pitman and Le (2005) to model lahar runout and dynamics. This model, implemented in the Titan2D toolkit (Pitman et al. 2003; Patra et al. 2005), simulates the flow of mixtures of soil, rock and interstitial fluid (water) over natural terrains. The main inputs required for this method is a digital elevation model (DEM), initial height and area of the mixture (pile height) and basal and internal friction values. The Pitman and Le (2005) model has previously been applied to simulate lahars (e.g. Williams et al. 2008; Procter et al. 2010; Córdoba et al. 2015). Crucially, this model does not consider the effects of erosion, deposition or transformation of the lahar and therefore will not reproduce all features of a lahar. Consequently, simulations using this method will not exactly match the outcomes of actual lahars, but will provide a probabilistic guide to the lahar hazard. In order to couple the susceptibility model material volumes to the inputs in the Titan2D toolkit, the Titan2D code was modified to add support for user defined pile heights in the form of a geographical information system (GIS) raster file. The modified source code is provided in the Zenodo repository https://doi.org/10.5281/zenodo.153993.

Study area and simulation inputs

The process outlined in the previous section is demonstrated through the generation of probabilistic hazard maps for lahars originating from the Mangatoetoenui Glacier, Ruapehu volcano, New Zealand. The Mangatoetoenui Glacier and Stream are located on the north-eastern flank of Ruapehu volcano, draining eastwards into Tongariro River (Fig. 1). A hydro-electric dam, trout fisheries and several towns located along Tongariro River are at risk from lahars originating on the north-eastern flanks of the volcano (Cronin et al. 1997). In this study, deposit characteristics are chosen to create a scenario with similar conditions to those in Mangatoetoenui Stream prior to the October 28, 1995 lahar described in Hodgson and Manville (1999) and used in lahar susceptibility estimations by Mead et al. (2016). The initial conditions in this study only consider a limited range of initial volumes and solid volume fractions, with all other inputs fixed, which means the generated hazard maps only represent a subset of the possible parameter space. A complete hazard assessment would consider the range of other inputs (such as basal and internal friction), their uncertainty and model suitability (Calder et al. 2015).

Fig. 1

Study area overview showing Ruapehu volcano (left), round the mountain track (purple), Mangatoetoenui Glacier (circled) and the North and South branches of Mangatoetoenui Stream (‘a’ and ‘b’ respectively)

Lahar simulations for were run for 30 min of simulation time using a 25 m resolution digital elevation model (DEM) sourced from Landcare Research NZ (2010). The DEM spanned the region from Mangatoetoenui Glacier (see Fig. 1) in the west to the confluence of Mangatoetoenui Stream and Tongariro River in the east, and was chosen because it best represented hydrological features of the terrain, despite the lower resolution compared to the contour-derived Land Information New Zealand DEM (Stevens et al. 2003).

Mobilised material volumes in the source study area were calculated for rainfall durations of 30 min, 2 h, 6 h, 12 h and 24 h at annual exceedance probabilities (AEP) of 0.5 (2 years annual recurrence interval (ARI)), 0.1 (10 years ARI), 0.02 (50 years ARI) and 0.01 (100 years ARI), using the susceptibility simulations first presented in Mead et al. (2016). We assumed most material is mobilised in a single event to match the observations of Hodgson and Manville (1999). GIS raster files of mobilised material depths from the susceptibility simulations were used as the pile height input in Titan2D. The internal friction angle of the material was given as 32°, which is between the minimum and maximum material angles of repose (Procter et al. 2010), basal friction was set to 20° to represent a reasonably fluid granular material.

In our simulations, the solids volume fraction of each modelled lahar was computed from the volumes of mobilised material and rainfall. The solid volume fraction was estimated as the volume fraction of solids in the total volume of mobilised material and water. Volume fractions, shown for each scenario in Table 1, decrease as ARI increases since higher rainfall intensities result in proportionally greater volumes of rainfall relative to the amount of mobilised material.

Table 1 Solid volume fractions of mobilised material for each scenario at Mangatoetoenui Stream

Generation of probabilistic lahar hazard maps

Lahar hazard maps derived from overall annual exceedance probabilities are shown in Fig. 2 for lahar height thresholds of 0.5 and 0.1 m. While lahar properties such as velocity and sediment concentration are also important to consider when evaluating overall lahar hazard, height is used here to express the hazard due to its ease in interpretation and use in delineation of hazard zones. Overall exceedance probabilities were calculated from simulated maximum lahar height at each grid cell using the complementary function

$$ 1-{p}_{\left( h\ngtr x\right)} $$

(1)

Fig. 2

Annual exceedance probabilities (AEP) for lahars exceeding 0.5 m and 0.1 m (bottom) in height

where p _{(h ≯ x)} is the probability of observing no lahar heights (h) greater than the threshold height x, defined using the formula of Tonini et al. (2016) as

$$ {p}_{\left( h\ngtr x\right)}={\prod}_i\left(1-{p}^i\left( h> x\right)\right) $$

(2)

where p ⁱ (h > x) is the probability of lahar heights exceeding the threshold for each AEP value i (in this instance i = 0.5, 0.1, 0.02, 0.01)

$$ {p}^i\left( h> x\right)= i\cdot p\left( h> x\right) $$

(3)

and p(h > x) is determined from the cumulative distribution of simulated lahar heights for each AEP.

Discussion and limitations

The deposit characteristics for this demonstration were chosen to create a scenario similar to the tephra deposit in the Mangatoetoenui Catchment prior to the October 28, 1995 lahars. A rainfall event (23 mm total in 24 h, 9 mm in one 6-h interval) triggered lahars in multiple catchments of Ruapehu (Manville et al. 2000), including the Mangatoetoenui catchment. The Mangatoetoenui lahar, described in Hodgson and Manville (1999), travelled downstream as a debris flow for the first 5 km (the proximal zone) and then progressively transformed into a hyperconcentrated flow between 5 and 9 km from the source due to the entrainment of streamflow and deposition of sediment. The lahar continued to dilute downstream, eventually reaching the Tongariro River. The lahar was confined to the Mangatoetoenui catchment, although a small deposit was observed in a distributary valley of the Whangaehu catchment (Hodgson and Manville 1999).

The similarities in initial conditions between the October 28 lahar and simulations mean that the hazard maps in Fig. 2 should show a reasonable degree of similarity between the hazard outlines and observations of the lahar. However; these hazard maps indicate that a large proportion of the lahar source material enters the Whangaehu catchment (shown in Fig. 2). Lahar heights also appear to fall below 0.1 m downstream of the walking track encircling Ruapehu (~8 km from source, purple line in Fig. 2). This differs from the extent of lahar deposits identified in Hodgson and Manville (1999), but is close to the distance (9 km) where flow fully transformed from a debris flow to hyperconcentrated flow. Differences in the DEM representation and initial deposit characteristics contribute to these observed errors, but the main source of error is the lack of representation of dilution, entrainment and flow transformation processes in current numerical models. The friction parameters used here provide a runout prediction near to the transition to hyperconcentrated flow, but also represent a highly mobile initial mass of material that causes a large proportion of the lahar volume to enter the Whangaehu catchment. This highlights the current state of lahar models as a key limitation of the proposed methodology and shows that operational hazard maps generated using this approach would still require expert guidance.

While differences in AEPs are visible in the vicinity of the walking track, there is little difference in lahar extent between AEPs in upstream portions of the Mangatoetoenui catchment (Fig. 2). This lack of difference can be attributed to the deeply incised channel walls present in upper reaches of Mangatoetoenui stream, but can also be caused by similar initial volumes of material being mobilised for all rainfall ARIs (see Mead et al. (2016) for explanation). However, the level of detail is limited in this demonstration as only 5 rainfall duration scenarios are simulated for each rainfall ARI, affecting the validity and resolution of the cumulative distribution used to determine event AEP. A complete hazard mapping exercise would need to sample the entire spectrum of rainfall durations in order to increase the resolution of the cumulative distribution.

Conclusion

The probabilistic lahar hazard maps generated here by combining numerical modelling with initial volumes determined through the susceptibility approach of Mead et al. (2016) have demonstrated a potential methodology for probabilistic hazard mapping of lahars. However, through this demonstration, some key limitations and simplifications have been identified that could affect the feasibility of this technique. Importantly, lahar numerical modelling approaches require further research to develop advanced methods capable of representing entrainment, deposition and flow transformations typical of lahars. Recent and ongoing research in this area (e.g. Pudasaini 2012; Pudasaini and Krautblatter 2014; Iverson and George 2016) may, in the future, provide alternative modelling approaches capable of accurately predicting additional features of lahar flows. Another limitation seen in this demonstration is the reduced parameter space used to generate the hazard maps. A complete lahar hazard assessment would need to quantify the range and uncertainty of all inputs into the model (e.g. basal and internal frictions) in addition to simulating a larger number of rainfall durations. However, simulations spanning the entire range of input parameters increases the computational requirements of this approach, which could impact on the feasibility of the method for lahar hazard assessment.