What drives landslide risk? Disaggregating risk analyses, an example from the Franz Josef Glacier and Fox Glacier valleys, New Zealand

We present a quantitative risk analysis (QRA) case study from the Kā Roimata o Hine Hukatere / Franz Josef Glacier and Te Moeka o Tuawe / Fox Glacier valleys, on the west coast of the South Island, Aotearoa / New Zealand. The glacier valleys are important tourist destinations that are subject to landslide hazards. Both valleys contain actively retreating glaciers; experience high rainfall; and are proximal to the Alpine Fault, which is a major source of seismic hazard on the west coast. We considered the life safety risk from rockfalls, soil/rock avalanches, and flows that either are seismically triggered or occur aseismically. To determine the range in risk values and dominant contributing variables to the risk, we modelled nine different risk scenarios where we incrementally changed the variables used in the risk model to account for the underlying uncertainty. The scenarios represent our central estimate of the risk, e.g. neither optimistic nor conservative, through to our upper estimate of the risk. We include in these estimates the impact time-variable factors, such as a recently reactivated landslide, have had on locally increasing risk and the time-elapsed since the last major earthquake on the nearby Alpine Fault. We disaggregated our risk results to determine the dominant drivers in landslide risk, which highlighted the importance of considering dynamic time-variable risk scenarios and the changing contributions to risk from aseismic versus seismic landslides. A detailed understanding of the drivers of landslide risk in each valley is important to determine the most efficient and appropriate risk management decisions.

1 Introduction

High-mountain areas are subject to a variety of natural hazards, including slope instability. Globally, these mountainous areas are currently experiencing declining low-elevation snow cover, retreating glaciers, and degrading permafrost as a result of climate change (see Hock et al., 2020). Such changes in environmental, meteorological, and geomorphological conditions may influence the rate, size, and characteristics of landslide hazards (Gariano and Guzzetti, 2016). Additionally, such high-relief mountain areas are subject to seismic hazards, including seismically triggered landslides. Given that the exposure of people and infrastructure to landslide hazards is also increasing from population growth, tourism, and socio-economic development (Hock et al., 2020), the risk from landslides may change and increase with time.

Quantitative risk analysis (QRA) is an important tool for assessing, managing, and communicating the risks from landslide hazards (Corominas et al., 2014), and there is an increasing need to undertake QRA coming from legislative authorities and from within the engineering and engineering geological communities (Corominas et al., 2014; Van Westen and Soeters, 2006; Ho et al., 2000). QRA is undertaken for land use planning (e.g. Bell and Glade, 2004; Vega and Hidalgo, 2016), infrastructure (e.g. Voumard et al., 2013; Macciotta et al., 2015), and visitor destinations (e.g. Corominas et al., 2019; Stock et al., 2014). However, the ability to estimate quantified levels of risk is often challenging as the input datasets used in risk analyses are inherently uncertain. Such uncertainty is mainly due to the lack of completeness, quality, or range within the input datasets required to undertake a QRA (Van Westen and Soeters, 2006). A landslide inventory, which details where landslides have occurred in the past, provides information critical for understanding what triggers landslides, what makes a particular slope more susceptible to landsliding, and how frequently landslides are likely to occur (Guzzetti et al., 2012). Yet, for many landslide-prone areas this spatial and temporal record of landsliding is limited or does not exist. This is particularly the case for certain trigger events, such as earthquakes, where the return period of the trigger event may be greater than the length of the historical record (van Westen et al., 2008). Consequently, assessments of landslide susceptibility and frequency rely heavily on practitioner experience and judgement (Lee, 2009) and may not always reveal the full levels of uncertainty attached to the risk estimates (Corominas et al., 2014; Macciotta et al., 2016). Most approaches use past landslide behaviour to predict what may occur in the future based on the maxim “the past is key to the future” (Varnes, 1978). However, the present or future conditions that make a slope susceptible to landsliding or trigger landsliding may be different to those of the past. Changes in the location and the frequency of landslide activity may substantially alter the estimated risk, adding to the uncertainty associated with the risk value.

Fell et al. (2005) suggest using sensitivity factor analysis as a tool to understand the influence of potential uncertainties on the estimated risk levels and to communicate the influence of this input variability to users of the risk analysis and assessment. A current limitation of risk analysis is the need to be able to “disaggregate” the risk results in order to determine the importance of the different input factors included in the QRA, such as the annual frequency of a given landslide type and volume occurring under a given set of triggers, how far landslide debris travels down a slope, where people are present on the slope, and their biophysical vulnerability if present and hit by landslide debris. Such limitation means that the contribution to the risk and sensitivity of the results relating to the input variables used is rarely quantified, thus making it difficult for risk managers to understand and implement targeted risk reduction measures and risk communication options. We address this here by presenting the QRA results and their uncertainties from local- to regional-scale (1:10 000–1:50 000) analyses of landslide hazards and the risk they pose to the lives of people visiting and working in the Kā Roimata o Hine Hukatere / Franz Josef Glacier and Te Moeka o Tuawe / Fox Glacier valleys, located on the west coast of Aotearoa / New Zealand. In the Southern Alps of New Zealand, landslides are a common feature that play a significant role in driving erosion (e.g. Hovius et al., 1997; Korup et al., 2004) and present an increasing natural hazard and risk to people and property (Allen et al., 2011; Cox et al., 2015; McSaveney, 2002).

Figure 1 Location of the Franz Josef Glacier and Fox Glacier valleys on the west coast of the South Island (a, b), which contains data sourced from the LINZ Data Service and licensed for reuse under the CC BY 4.0 licence. (c) Photograph of the Franz Josef Glacier valley. (d) The Fox Glacier valley, including the access roads, cycleway, and tracks within it. (e) The Franz Josef Glacier valley, including the access roads, cycleway, and tracks within it.

2 Study site

The glacier valleys, which are important tourist destinations, are located on the west coast of the South Island, New Zealand (Fig. 1). Both valleys contain multiple trails (walking and/or cycling tracks), which take between 30 min and 8 h to walk/cycle, that allow visitors to access and experience a glacial environment. The glaciers themselves can now only be accessed via helicopters, with visitors undertaking paid tours on the glacier. Given the commercial sensitivities, the risk from landslide hazards to visitors on commercial tours on the glacier has not been quantified in our risk analysis. Prior to the Covid-19 pandemic and associated closure of New Zealand's border to international tourists, ca. 700 000 people per year walked the tracks in the Franz Josef Glacier valley and ca. 400 000 people per year walked the tracks in the Fox Glacier valley. A maximum of 6000 and 3500 people per day walked the tracks in the Franz Josef Glacier valley and Fox Glacier valley, respectively. Within this environment, visitors are exposed to a variety of landslide hazards. Numerous near misses have been documented, and two fatalities occurred in January 1980 when a debris avalanche occurred along a track in the Fox Glacier valley. Currently, the northern road and access track within the Fox Glacier valley are closed due to repeated damage from debris flow events. Evidence of landsliding is present within each valley, with the types of landslide broadly classified into rockfalls, slides and topples, debris and rock avalanches, and debris flows (as classified by Hungr et al., 2014; see Fig. 2). In addition to these broad landslide types, deep-seated gravitational slope deformations (DSGSDs) can be observed in both study areas. These large DSGSDs typically provide sources of material for smaller rockfall/debris avalanches or debris flows (Cody et al., 2020). Earthquakes are potential triggering mechanisms for landslides as both study areas are located less than 10 km south-east of the Alpine Fault (Fig. 1), which is a major source of earthquakes in New Zealand (Stirling et al., 2012). Additionally, both valleys experience high rainfall, with 5 m yr⁻¹ recorded in the village of Waiau / Franz Josef increasing to >10 m yr⁻¹ towards the main divide of the Southern Alps (Langridge et al., 2016). The glaciers in each valley are currently retreating (Purdie et al., 2015, 2021), exposing more disturbed and consequently weaker rock masses, which appear to be the source of many recently documented landslides. Many of these “aseismic” landslides appear to be triggered by intense rainfall; however, several have no documented trigger. Therefore, the slopes in the study areas have been and will continue to be subjected to transient changes in stress, typically caused by precipitation-induced variations in pore-water pressure, erosion, freeze–thaw cycles, and diurnal and seasonal temperature variations (Gunzburger et al., 2005; Viles, 2013; Eppes and Keanini, 2017). Transient stress changes within a slope can lead to deformation, fracturing, and joint dilation, thus reducing rock mass strength and leading to failure (e.g. Eberhardt et al., 2004; Eppes and Keanini, 2017).

The study areas (shown in Fig. 1) are dominated by ice-free slopes comprised of schist (Cox and Barrell, 2007). The structural geology of the bedrock schist is complex, given the proximity of both sites to the Alpine Fault (Fig. 1). Large persistent faults cut through the area trending north-east to south-west and east to west. The quality of the rock mass is highly variable over the study areas and tends to change with proximity and location relative to these persistent faults. Moraine and colluvium deposits are present within the main and tributary valleys, with the valley floors formed of predominantly alluvium and re-worked moraine and colluvium. The glaciers have carved the valleys, resulting in steep bedrock valley sides truncated by deeply incised streams. Debris fans are present at the mouth of these incised streams, which feed into the Fox and Waiho (Franz Josef) rivers. Flow rates within these rivers are highly seasonal, and their courses, within their respective wider valleys, change frequently.

In the Fox Glacier valley, there are more extensive and thicker debris deposits (both moraine and colluvium) and larger debris fan deposits (e.g. Yellow Creek fan; Gomez and Purdie, 2018), indicating that debris flows and avalanches may be more prevalent. With glacier retreat, these debris deposits are free to begin creeping and debris is available for remobilisation via debris flows (Cody et al., 2020). In contrast, Franz Josef contains less debris and is more dominated by bedrock slopes, which may be the result of limited debris accumulation through time or the ability of erosional processes to keep pace with deposition and remove material from the valley.

3 Methodology

3.1 Risk calculation route

To estimate the risk, we follow the quantified risk analyses method described in AGS (2007) and the Joint Technical Committee on Natural Slopes and Landslides (JTC1, as outlined in Fell et al., 2008). We calculated the probability of death (life risk) of an individual, P_(LOL), from

$$\begin{array}{}\text{(1)}& P_{\left(\mathrm{LOL}\right)}=P_{\left(\mathrm{L}\right)}\times P_{(\mathrm{T}:\mathrm{L})}\times P_{(\mathrm{S}:\mathrm{T})}\times V_{(\mathrm{D}:\mathrm{T})},\end{array}$$

where P_(L) is the probability (annual frequency) of the landslide occurring. P_(T:L) is the probability of the landslide (e.g. the debris from a landslide of a given type) reaching the element at risk (e.g. visitor on a track). P_(S:T) is the spatio-temporal probability of the person at risk being present and in the path of landslides (the proportion of a year that the person is exposed to landslides). V_(D:T) is the vulnerability of the person if present and in the path of landslide debris (i.e. the probability that the person will be killed if impacted by the landslide). In our analyses, we include in the vulnerability estimates the potential for a person to be aware of the hazard and take evasive action.

Within the Fox Glacier and Franz Josef Glacier valleys, we sum risk from several landslide hazards where people are exposed to (1) several types of landslides, (2) landslides of the same type but different volume, (3) landslides triggered by more than one phenomenon, and (4) several slopes on which landslides can occur. To take such cases into account, we re-wrote Eq. (1) as

$$\begin{array}{}\text{(2)}& P_{\left(\mathrm{LOL}\right)}=\sum _{i=\mathrm{1}}^n\left(P_{i\left(\mathrm{L}\right)}\times P_{i(\mathrm{T}:\mathrm{L})}\times P_{i(\mathrm{S}:\mathrm{T})}\right)\times V_{i(\mathrm{D}:\mathrm{T})},\end{array}$$

where n is the number of landslide hazards of a given type and volume. This assumes that the hazards are independent of each other. However, in the valleys, it is possible that one or more of the hazards may result from the same causative event, e.g. an earthquake. Therefore, we estimate the probabilities using the theory of uni-modal bounds where the upper-bound conditional probability (P_UB) is calculated from

$$\begin{array}{}\text{(3)}& P_{\mathrm{UB}}=\mathrm{1}-\left(\mathrm{1}-P_{\mathrm{1}}\right)\times \left(\mathrm{1}-P_{\mathrm{2}}\right),\mathrm{\dots },\left(\mathrm{1}-P_n\right),\end{array}$$

where P₁ to P_n are the estimate of several individual hazard conditional probabilities. We then multiplied the P_UB by the annual probability of the common causative event, e.g. the given level of shaking representing a given earthquake. More detail on the equation route is provided in the Appendix.

3.2 Risk metrics

We estimated the risk to life using four risk metrics. Firstly, we estimated the local personal risk (LPR), which represents the annual probability of death for a hypothetical person present at a particular location for 100 % of the time (24 h d⁻¹ and 365 d of the year). LPR, a metric used in flooding and seismic hazard studies (Crowley, 2017; Jonkman et al., 2003; van Elk et al., 2019), can be used to visualise the spatial distribution of risk within the study areas in order to help plan/realign tracks and roads. Secondly, we estimated the individual risk per trip, which is expressed in terms of the fatality risk (probability of death) of an individual resulting from one return trip along one of the main access tracks or roads within the study areas. We use this to represent the risk to visitors. We then estimated the annual individual fatality risk (AIFR), which is expressed in terms of the fatality risk experienced by the most exposed individual over 1 full year of, for example, working in the valleys and undertaking frequent track checks for 2 to 3 h d⁻¹ (see Massey et al., 2022a). We use AIFR to estimate the risk to the most exposed worker in each valley who is present every day for substantial periods of the year. We estimated societal risk by determining f–N pairs, which represents the frequency (f) of an accident killing one or more people (N) in a single event, plotted on an f–N curve (Strouth and Mcdougall, 2021). In this paper, we focus on and report the results for LPR and individual risk per trip. AIFR and the f–N curve results are reported in Massey et al. (2022a).

3.3 Methodology framework

Our risk analysis firstly considers the possible range of triggering events in terms of a set (bands) of earthquake triggers and aseismic triggers (e.g. rain, time). In our compilation of the landslide inventories for each valley, we were unable to determine a relationship between rainfall or snowmelt with landslide occurrence. The recorded near misses in the Franz Josef Glacier valley and two fatalities (January 1980) in the Fox Glacier valley from a debris avalanche occurred in the absence of any discernible trigger. Therefore, we subsume potential rainfall triggering, snowmelt triggering, and climatic factors into an aseismic annualised rate of landsliding. However, due to the proximity of the Alpine Fault and the seismic history of the region, we explicitly considered the possibility of seismically triggered landslides.

Figure 2 Schematic and example photographs of landslide types considered in the QRA, including (a) rockfall, (b) debris and rock avalanches, and (c) debris flows.

For each representative earthquake event, we determined the annual frequency of the event and the number of landslides of a given volume class produced in that event. For aseismic landslides, we determined the annual frequency of landslides of a given volume occurring in each valley using historical data on aseismic landslides in the valleys and the wider Southern Alps (Massey et al., 2022a). For both seismic and aseismic landslides we considered the full range of volume classes that could occur in each valley, which are (1) ≤10 m³, (2) 10 to 100 m³, (3) 100 to 1000 m³, (4) 1000 to 10⁴ m³, (5) 10⁴ to 5×10⁴ m³, (6) 5×10⁴ to 10⁵ m³, (7) 10⁵ to 5×10⁵ m³, (8) 5×10⁵ to 10⁶ m³, (9) 10⁶ to 5×10⁶ m³, and (10) $>\mathrm{5}\times {\mathrm{10}}^{\mathrm{6}}$ m³. We estimated the number of landslides that could occur for each volume class using the Moon et al. (2005) method, by calculating the area under the landslide volume–frequency curve (see Fig. 3) using log–log histogram bins.

Secondly, we considered the locations from which landslides are most likely to be sourced in each glacier valley. We explicitly determined landslide source locations in order to estimate how far the debris could travel downslope from a particular source. We used slope angles, volume-to-area scaling relationships, and geomorphic mapping to delineate these source areas. We compiled information on pre-disposing factors of slope instability in each valley to understand spatial controls on landslide occurrence, with these datasets forming an important input into landslide susceptibility modelling for each valley (Reichenbach et al., 2018). We used logistic regression susceptibility models for both seismic (Massey et al., 2021) and aseismic landslides to weight which source areas may preferentially generate landslides.

We conducted 3D numerical runout simulations to determine P_(T:L): the probability of the debris from a landslide reaching or passing a portion of slope as it travels downhill from the source area. We conducted these numerical simulations for rockfall, debris avalanches, and debris flows from our explicitly determined source areas. We used RAMMS rockfall software (RAMMS, 2015) for rockfall simulations and RAMMS debris flow software (RAMMS, 2011) for debris flow and debris avalanche simulations.

We compiled information on the length of time visitors and workers spend along the tracks in each valley to estimate the spatio-temporal probability of the person at risk being present at a location (P_(S:T)) and consequently in the path of debris. We used empirical estimates of vulnerability (V), which is the probability of a person being killed if present and in the path of one or more boulders, considering both (a) the likelihood of being killed if struck and (b) the possibility of being able to take evasive action and avoid being struck.

For each of these steps and elements of the risk equation, we determined central estimates (which we define as neither optimistic nor conservative and are based on taking the mean estimate) and upper estimates (based on the 84th percentile) of the different variables used in the QRA. We used these different estimates in our sensitivity analysis, where we incrementally changed the variable estimates from central to upper until all variables were upper estimates. This results in nine different risk models, from which we can calculate the incremental change in risk based on varying the input assumptions and document the impact of aleatory uncertainty.

3.4 The probability (annual frequency) of the landslide occurring – P_(L)

3.4.1 Seismic landslides

We determined the frequency and volume of landslides likely to be generated at different magnitudes of ground-shaking intensity from the mapped landslide distributions of historical New Zealand and international earthquakes, as detailed in de Vilder et al. (2020). We used the landslides generated during the 2016 M_w 7.8 Kaikōura, 1968 M_w 7.1 Inangahua and 1929 M_w 7.8 Murchison earthquakes as proxies (Massey et al., 2018; Hancox et al., 2014, 2015). We selected these three landslide inventories as they represent the most complete New Zealand inventories for seismic landslides that occurred in fractured hard rock (such as greywacke) similar to that of schist and occurred in mountainous and hilly terrain.

We assessed the number of landslides that could be generated for four different representative earthquake events, as represented by peak ground acceleration (PGA) bands: Band 1 (0.2–0.35 g), Band 2 (0.35–0.65 g), Band 3 (0.65–1.2 g), and Band 4 (≥1.2 g). It is unlikely that several landslides would be generated by ground shaking < 0.2 g (Dowrick et al., 2008). We calculated the annual frequency of the representative PGA per band from the New Zealand National Seismic Hazard Model (NSHM) (Stirling et al., 2012) by subtracting the annual frequencies that represent the PGA boundaries (start and end) of each band. The active fault component of the NSHM defines the Alpine Fault local to the Franz Josef Glacier valley and Fox Glacier valley as the AlpineF2K fault source. Within the NSHM the AlpineF2K source generates a M_w 8.1±0.2 earthquake with a single-event (strike-slip + dip-slip) displacement of ca. 9.2 m with a mean recurrence interval of 341 years (Stirling et al., 2012). This is a time-independent variable and does not consider the time elapsed since the last earthquake on the Alpine Fault in 1717 (Howarth et al., 2021). Landgride et al. (2016) disaggregated the NSHM to see what other fault sources may contribute to the shaking hazard at Franz Josef. For a probability of roughly 10 % in 250 years (or 2500 years) the disaggregation indicates that the main contributor of seismic hazard is the M_w 8.1 AlpineF2K source (i.e. the Alpine Fault). Additionally, the second-largest seismic hazard over 2500 years comes from moderate-magnitude (M_w 5–6) earthquakes that can occur <10 km from the townships. Although the Alpine Fault is the main seismic source in the area, the section of fault that could rupture might be located some distance away from the sites. For this reason and to consider the contribution from the M_w 5–6 earthquakes, we therefore estimate the landslide severity for the four different bands of PGA as determined from the NHSM. Recent research (Howarth et al., 2021) shows that the probability of an earthquake occurring on the central section of the Alpine Fault is 75 % in the next 50 years and that there is an 82 % chance that the earthquake will be greater than M_w 8. To account for this in our upper estimate of the PGA annual frequencies, we increase the annual frequency of the most intense ground shaking (Band 4) to 0.015 to reflect the time elapsed since the last Alpine Fault earthquake.

To assess the magnitude–frequency of seismic landslides in each band, as outlined in de Vilder et al. (2020), we firstly determined the appropriate landslide source volume-to-area scaling relationship (from Massey et al., 2020). Secondly, we estimated relationship between the landslide frequency (number) and source area scaling. Thirdly, we investigated the relationship between landslide occurrence and PGA, slope angle, and material type using the Kaikōura, Inangahua, and Murchison landslide inventories. Finally, we combined estimates of the annual frequency of the representative event PGA for each earthquake band in the NSHM. Using this relationship, we estimated the probability of a landslide of a given volume class occurring within each study area for each PGA band considered, along with the annual frequency of the representative PGA in the band occurring (see Fig. 3a). We fitted power laws to the data, with these representing our central estimate of the number of landslides of a given volume class occurring for each PGA band considered. To derive an upper estimate, we added the standard error of the gradient of our best-fit power law to the power-law relationship to calculate the number of landslides that could be generated.

Figure 3 (a) The number of landslides of given landslide volume that could be generated for different levels of ground shaking as represented by Bands 1 through 4. The power law for Band 2 is $y=\mathrm{2.28}\times {\mathrm{10}}^{\mathrm{6}}\times \left({\mathrm{vol}}^{-\mathrm{1.60}}\right)$, for Band 3 is $y=\mathrm{2.33}\times {\mathrm{10}}^{\mathrm{6}}\times \left({\mathrm{vol}}^{-\mathrm{1.45}}\right)$, and for Band 4 is $y=\mathrm{1.47}\times {\mathrm{10}}^{\mathrm{7}}\times \left({\mathrm{vol}}^{-\mathrm{1.41}}\right)$. (b) Magnitude–frequency relationships for aseismic landslides for the Fox Glacier valley and Franz Josef Glacier valley. The two power laws on each graph represent the central estimate (using the power-law relationships). The power law for the Franz Josef Glacier valley is $y=\mathrm{9.08}\times \left({\mathrm{vol}}^{-\mathrm{0.66}}\right)$, and for the Fox Glacier valley it is $y=\mathrm{5.21}\times \left({\mathrm{vol}}^{-\mathrm{0.51}}\right)$.

3.4.2 Aseismic landslides

We collated information on the occurrence of aseismic landslides from various data sources, which we used to assess the historical type, mechanisms, and rates of aseismic landslides for both valleys. These data sources include (1) a rockfall register compiled from observations made by staff of the New Zealand Department of Conservation (DOC), Franz Josef Glacier Guides Ltd, and Fox Glacier Guiding; (2) a landslide inventory derived from historical aerial imagery analysis of both valleys; and (3) a large landslide inventory of historical landslides observed in the wider Southern Alps (see the Appendix for more information on the compilation of the landslide inventories).

We determined valley-specific magnitude–frequency relationships of landslides, given the amount of catchment-specific information about landslides. We fitted power-law trends to the data to generate a central estimate and an upper estimate (using the standard error of the power law) of the number of landslides that could occur in each valley per square kilometre and their annual frequency (see Fig. 3b). These landslide rates were then scaled to each valley by multiplying them with the total area of the slopes greater than 30^∘ within each valley. We used a slope angle of 30^∘ as a cut-off as within our landslide inventory no landslides occurred on slopes with angles < 30^∘.

3.5 The probability of the landslide reaching the track or road

3.5.1 Landslide susceptibility

In the absence of historical information on landslides triggered by an Alpine Fault earthquake, we used the 2016 M_w 7.8 Kaikōura earthquake-induced landslide inventory to understand the spatial controls on susceptibility to failure. From this dataset, Massey et al. (2018, 2021) used a logistic regression model to correlate the three aforementioned mapped earthquake-induced landslide inventories with various topographic, geological, and seismological parameters to understand which parameters best explained the occurrence of coseismic landslides. We applied the Massey et al. (2018) logistic regression model to both valleys, using PGA input from the NSHM. The PGA input varied from 0.8 to 1.1 g across both valleys, which, along with the other components of the Massey et al. (2018) logistic regression model such as distance to fault, slope angle, geology, and local slope relief, was used to determine the probability of a landslide occurring from a particular source location.

We developed valley-specific logistic regression models to determine aseismic landslide susceptibility, given the amount of landslide-specific information and slight differences in the landslide hazards within each valley (see Fig. 4 for an example from the Fox Glacier valley). These models are based on the correlation of mapped landslides and various topographic, geological and land use characteristics (see Massey et al., 2018, and the Appendix). Rockfalls recorded in the rockfall register were not included within the analysis as the data do not have accurate geographic locations. More information on the aseismic susceptibility models is provided in the Appendix.

Figure 4 Generation of both aseismic and seismic logistic regression models for the Fox Glacier valley. The aseismic logistic regression model is determined from the correlation between (a) the landslide inventory and the static variables of (b) slope angle, (c) local slope relief (LSR, defined as the maximum height difference within a fixed 80 m radius of the grid cell), (d) geology, and (e) vegetation (as classified from imagery and the New Zealand Land Cover Database). The logistic regression model calculates (f) aseismic landslide susceptibility. The seismic logistic regression model is determined from the correlation between (g) peak ground acceleration (PGA) and the static variables of (b) slope angle, (c) local slope relief (LSR), (d) geology, and (h) distance to active faults (Massey et al., 2018). The logistic regression model calculates (i) seismic landslide susceptibility.

3.5.2 Landslide runout

Landslide sources ≤ 1000 m³ were assumed to be rockfalls rather than debris flows and debris avalanches. Potential rockfall source areas were defined using all slopes ≥ 45^∘, assuming any slope ≥ 45^∘ can potentially generate rockfalls (Fig. 5) (Budetta, 2010; Massey et al., 2014a). For landslide volumes ≤ 10⁵ m³ the landslide sources were assumed to be pixels of a given area based on the area-to-volume scaling exponents (Fig. 5). For landslide volumes >10⁵ m³, the shapes of the sources were defined using the geomorphic features (Fig. 5).

Figure 5 (a) Simulated source areas in the Fox Glacier valley for all volume classes. (b) Numerical rockfall simulations, using RAMMS, from pixel source areas. (c) Numerical debris avalanche and debris flow simulations, using RAMMS, from pixel source areas for volume classes 10⁴ to 10⁵ m³. (d) Numerical debris avalanche simulations, using RAMMS, for geomorphically defined source areas for volume classes $\ge \mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ m³.

For the rockfall simulations, we used RAMMS rockfall software (RAMMS, 2015), which simulates the rigid-body motion of falling rocks and predicts rock trajectories in general three-dimensional terrain (see the Appendix for more information). For the debris avalanches and debris flow simulations, we used RAMMS debris flow software (RAMMS, 2011), changing the Voellmy friction parameters (see the Appendix for more information) to determine if a particular source area failed as a debris avalanche or debris flow. From these simulations, we derived the runout extent and maximum debris height. The simulated maximum height of debris passing through a given grid cell is converted into the number of boulders (our central estimate), with our field measurements indicating that the median boulder size is 1 m³ (Fig. 6). For example, if the maximum debris height passing through a 3 m by 3 m grid cell is 1 m, then the total volume of debris passing through that grid cell is 9 m³, which when converted into N boulders, would be on average nine boulders. Based on sensitivity analysis of the Voellmy friction parameters (see the Appendix for more information), we calculated a standard-deviation-based factor of difference in debris height. We applied this factor of difference to the simulation results to increase debris height, providing an upper estimate.

Figure 6 (a) Boxplot of the measured boulder volumes (n=36) in the field, indicating a median boulder volume of 1 m³. (b) Histogram of the measured boulder volumes in the field.

We calculated the probability of one boulder in the debris hitting an object when passing through a particular portion of the slope, perpendicular to the debris path, using the following equation:

$$\begin{array}{}\text{(4)}& P_{\mathrm{1}(\mathrm{T}:\mathrm{L})}={\displaystyle \frac{D+d}{L}},\end{array}$$

where D is the diameter of the boulder, d is the diameter of a person (our central estimate assumes a person is a “cylinder” with a 1 m diameter, while our upper estimate assumes a 2 m diameter), and L is the unit length of slope perpendicular to the runout path, which for this study is 3 m grid cell. Our equation includes a “buffer zone” around the person (D+d) within which the boulder travels along a path either side of d and cannot miss.

3.6 Exposure

We calculated the probability that a person will be occupying a given grid cell along one of the tracks/roads (P_(S:T)) if they spend a number of hours (N_HRS) per trip per year walking/driving that route using the following equation:

$$\begin{array}{}\text{(5)}& P_{(\mathrm{S}:\mathrm{T})}={\displaystyle \frac{\left(N_{\mathrm{HRS}}\right)}{\left(N_{\mathrm{C}}\right)}},\end{array}$$

where N_C is the number of cells visited along the route. We compiled information provided by DOC on the estimated time taken to travel by vehicle along given roads to/from the car parks and the time taken for walking a round trip from the car parks to the glacier viewing points to determine the time exposed for an average walker (central estimate) or a slow walker (upper estimate). In the Fox Glacier valley, our average walker spent 1.5 h walking to and from the glacier viewpoint and 0.2 h driving to and from the car park (see Fig. 1c), while the slower walker spent 2 h walking to and from the glacier viewpoint and 0.3 h driving to and from the car park. In the Franz Josef Glacier valley, our average walker spent 2 h walking to and from the glacier viewpoint and 0.3 h driving to and from the car park (see Fig. 1d), while the slower walker spent 2.5 h walking to and from the glacier viewpoint and 0.4 h driving to and from the car park. Therefore, we assumed that the time spent on each 1 m section of track was equal to the duration (time) of travel divided by the total length of the track. However, exposure can be adjusted to account for a longer time spent by a visitor at a viewing area, picnic spots, etc. For the calculation of LPR, we assumed a P_(S:T) of 1, where a person is present 100 % of the time.

3.7 Vulnerability

Physical vulnerability (V) depends on the landslide intensity, the characteristics of the elements at risk, and the impact of the landslide (Du et al., 2013). To derive our central estimate of vulnerability, we link vulnerability values to representative landslide volumes, which act as a proxy for landslide intensity (see the Appendix for more information). Anecdotal evidence from the glacier valleys suggests that evasive action reduces vulnerability. This was the case on 13 October 2011 and 16 June 2014, when boulders and flyrock from debris avalanches passed over several people on the ice. During these near misses, the guides heard the debris moving down the slope and had time to instruct their clients to take evasive action. In Table 1 our vulnerability values scale with landslide volume, where for landslide volumes ≤10⁵ m³ an individual may be able to take evasive action. However, the ability to take evasive action decreases with landslide volume. For landslide volumes > 10⁵ m³ an individual is likely to be buried by debris and killed. To derive an upper estimate of vulnerability, we assume a vulnerability of 1 for all rockfall and landslide volumes.

Table 1 Physical vulnerability values used in study.

3.8 Rick scenarios modelled

For each valley, we estimated the individual risk per trip for a visitor using nine different risk model scenarios (Table 2) which ranged from our central estimate of the risk to our upper estimate of the risk. Our central-estimate (Scenario 1 – Table 2) risk model uses the central-estimate input variables and a time-independent ground-shaking annual frequency. For each risk model scenario we incrementally change the variable from central to upper estimates, starting with the number of landslides that could occur in an earthquake event (Scenario 2) and then the number of aseismic landslides that occur annually (Scenario 3), before increasing our estimate of debris height (Scenario 4), the diameter of a person (Scenario 5), our vulnerability estimate (Scenario 6), the time elapsed since the last Alpine Fault earthquake (Scenario 7), and lastly the length of time a visitor is exposed to the risk (Scenario 8). We also included a risk model scenario (Scenario 9 – Table 2), where we used central-estimate input variables but account for the increased probability of Alpine Fault earthquake occurring to understand the impact of these assumptions on the risk results.

Table 2 Risk model scenarios and associated input variables.

4 Results

4.1 Individual risk per trip

The individual risk per trip in the Franz Josef Glacier valley ranged from $\mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{7}}$ (central estimate – Scenario 1) to $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{6}}$ (upper estimate – Scenario 8). The risk along the road in Franz Josef ranges from $\mathrm{7.88}\times {\mathrm{10}}^{-\mathrm{9}}$ to $\mathrm{1.03}\times {\mathrm{10}}^{-\mathrm{7}}$, while the risk along the track in Franz Josef ranges from $\mathrm{7.72}\times {\mathrm{10}}^{-\mathrm{7}}$ to $\mathrm{1.08}\times {\mathrm{10}}^{-\mathrm{5}}$. For Scenario 9 (central estimate with higher earthquake annual frequency), the risk along the road in Franz Josef was $\mathrm{6.84}\times {\mathrm{10}}^{-\mathrm{8}}$ and the risk along the track was $\mathrm{2.73}\times {\mathrm{10}}^{-\mathrm{6}}$, with a total risk per trip of $\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{6}}$. The individual risk per trip in the Fox Glacier valley ranged from $\mathrm{4.9}\times {\mathrm{10}}^{-\mathrm{6}}$ (central estimate – Scenario 1) to $\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{5}}$ (upper estimate – Scenario 8). The risk along the road in Fox ranged from $\mathrm{2.57}\times {\mathrm{10}}^{-\mathrm{7}}$ to $\mathrm{6.34}\times {\mathrm{10}}^{-\mathrm{7}}$, while the risk along the track in Fox ranged from $\mathrm{4.63}\times {\mathrm{10}}^{-\mathrm{6}}$ to $\mathrm{1.62}\times {\mathrm{10}}^{-\mathrm{5}}$. For Scenario 9 (central estimate with higher earthquake annual frequency), the risk along the road in Fox was $\mathrm{4.22}\times {\mathrm{10}}^{-\mathrm{7}}$ and the risk along the track was $\mathrm{7.16}\times {\mathrm{10}}^{-\mathrm{6}}$, with a total risk per trip of $\mathrm{7.59}\times {\mathrm{10}}^{-\mathrm{6}}$. The risk along roads is less than that along tracks; this is a function of both overall lower LPR risk and less time spent on the roads. It is important to note that the risk numbers reported here do not consider any risk management and mitigation so should not be treated as indicative of current residual risk levels following actions taken in light of this analysis.

Figure 7 Risk results from Scenario 1 (central estimate) risk model for the Franz Josef Glacier valley. (a) LPR map displaying areas of higher and lower risk, along with the location of tracks (black lines) in the valley. The risk results presented in (b) to (d) are extracted along the lower track in (a). (b) LPR values for aseismic landslides compared with seismic landslides. (c) LPR values for different volume classes of aseismic landslides. (d) LPR values for the different bands of ground shaking (from the lowest, Band 1, to the highest, Band 4).

Figure 8 Risk results from Scenario 1 (central estimate) risk model for the Fox Glacier valley. (a) LPR map displaying areas of higher and lower risk, along with the location of the access track (black line) in the valley. The risk results presented in (b) to (d) are extracted along the track in (a). (b) LPR values for aseismic landslides compared with seismic landslides. (c) LPR values for different volume classes of aseismic landslides. (d) LPR values for the different bands of ground shaking (from the lowest, Band 1, through to the highest, Band 4).

Figure 9 Risk results for Scenario 9 (central estimate and time-dependent earthquake frequency) along the access track in the (a) Franz Josef Glacier valley and (b) Fox Glacier valley, displaying the contributions of aseismic landslide risk and seismic landslide risk.

4.2 Risk disaggregation

Using Scenario 1, we disaggregate our risk results to understand the contributions to risk from the different risk model components. Figures 7 and 8 display an LPR map of Franz Josef and Fox, respectively, illustrating the spatial variation in risk within the valley and along the access tracks to the viewpoint of the glaciers. Aseismic landslides account for 66 % and 83 % of total LPR along the access tracks in Franz Josef and Fox, respectively, in contrast to 34 % and 17 % for seismic landslides (Figs. 7b and 8b). In Fox, increases in aseismic landslide risk are observed when the track is close to the base of larger, steeper slopes or crosses a large debris fan (Fig. 8a). Although the risk to an individual is higher for aseismic landslides, the risk of a large landslide causing multiple fatalities or multiple landslides occurring at the same time leading to multiple fatalities is dominated by earthquake events with PGA values > 0.6 g (Band 3 and Band 4) (see Massey et al., 2022a). For aseismic landslides, we disaggregated the risk further to determine which landslide volume classes contributed most to the risk. In Franz Josef, moderately sized landslide volume classes of 10⁴, 5×10⁴, and 10⁵ m³ account for 31 %, 30 %, and 14 % of the aseismic landslide risk along the track, while landslide volume classes of 5×10⁵ and 10⁶ m³ account for a further 7 % and 11 % of LPR, respectively (Fig. 7c). Landslide volume classes of 5×10⁶ m³ or greater account for less than 3 % of LPR. The risk from rockfalls (4 % of LPR along track) increases when the track is closer to the base of the steep valley sides as displayed in the spikes in risk associated with volume classes of 10 m³ (Fig. 7c). Increases in the risk associated with the 10⁴ m³ landslide volume is associated with increases in both aseismic and seismic risk along the track (Fig. 7). In the Fox Glacier valley, landslide volume classes of 10⁴, 5×10⁴, and 10⁵ m³ account for 24 %, 21 %, and 12 % of LPR, respectively (Fig. 8c). Larger volume classes of 5×10⁵, 10⁶, 5×10⁶, and $>\mathrm{5}\times {\mathrm{10}}^{\mathrm{6}}$ m³ contribute 12 %, 15 %, 10 %, and 6 % to LPR, respectively. For seismic landslides in Franz Josef, Band 2 contributes the most to the risk, accounting for 43 % of LPR, while Band 1 accounts for 21 %, Band 3 for 32 %, and Band 4 for 4 % of LPR (Fig. 7d). A similar pattern exists for seismic landslides in Fox, where Band 2 contributes the most risk, accounting for 48 % of LPR, while Band 1 accounts for 29 %, Band 3 accounts for 21 %, and Band 4 accounts for 2 % of LPR (Fig. 8d).

However, in Scenario 9 (Table 2), we modelled the increased annual frequency of a large Alpine Fault event that was assumed to result in the greatest ground shaking (Band 4) while using the central estimate for all other input variables to the risk model. In this scenario, the seismic landslide risk in Franz Josef (Fig. 9a) is higher, accounting for 81 % of LPR along the track, than that of aseismic landslide risk, in contrast to the patterns in Scenario 1 (Fig. 7b). In Fox, the contribution of seismic landslide risk is higher in Scenario 9, accounting for 46 % of LPR along the track (Fig. 9b), than in Scenario 1 (Fig. 8b), and in locations along the track, it surpasses that of aseismic landslides. However, aseismic landslides contribute more to the overall LPR (54 %), particularly in locations where the track crosses debris fans (Fig. 9b).

Figure 10 (a, b) Photographs of the evolution of the Mill's Creek debris fan and Alpine Gardens landslide. (c) The linkage between the Alpine Gardens landslide and the Mill's Creek debris channel and fan.

For eight different risk scenarios (Table 3), we calculate the overall cumulative increase in risk as a percentage, along with the amount of cumulative increase in risk between each scenario. We also calculate the increase in risk between each scenario as a percentage to understand the contribution of each variable to overall risk and negate the effect of the order in which each variable is altered in the risk models and the compounding effect of changes to variables on the cumulative risk results. Our sensitivity analysis of the risk model inputs for Franz Josef (Table 3) shows that, within the ranges of inputs considered, the largest increase in the risk is associated with increasing the number of aseismic landslides that occur annually with a cumulative increase in risk of 365 %. Second to this is the increased earthquake annual frequency, with a cumulative increase in risk of 330 %. Thirdly, increased exposure time results in a cumulative risk increase of 328 %, while a constant vulnerability of 1 results in a cumulative increase in risk of 219 %. Changes in input variables (debris height and diameter of a person) that affect the P_(T:L) term resulted in negligible changes to risk, with <10 % change in cumulative risk. Overall, changes in the input variables from the central estimate to upper estimate resulted in a 1298 % cumulative increase (just over an order of magnitude) in the risk results. For Fox, increases in exposure time and vulnerability resulted in the largest increase in risk (increase in cumulative risk of 80 % and 60 %, respectively; Table 3), while changes in the annual frequency of earthquake events and the number of aseismic landslides resulted in cumulative increases in risk of 56 % and 30 %, respectively (Table 3). Changes in the number of seismic landslides resulted in a cumulative increase in risk of 10 %. Similarly to Franz Josef, changes in debris height and the diameter of a person had negligible impact (Table 3). The range in risk values for Fox from the central estimate to upper estimate was smaller than for Franz Josef, with a cumulative percentage increase of 244 %.

Table 3 Risk model sensitivity analysis displaying the factor in increased risk between each scenario and the cumulative increase in risk from the central to upper estimate (Scenario 8).

n/a: not applicable.

4.3 Changing risk through time

Recently within the Fox Glacier valley, there has been increased debris flow events from the Mill's Creek catchment. The debris flows are sourced from the toe of the Alpine Gardens landslide, which is an approximately $\sim \mathrm{50}\times {\mathrm{10}}^{\mathrm{6}}$ m³ actively moving landslide complex in the Fox Glacier valley (Fig. 10). These debris flows travel down Mill's Creek and deposit on the debris fan at its confluence with the Fox River. The debris flow activity has resulted in the expansion of the Mill's Creek debris fan, which, in turn, has forced the migration of the Fox River to the true right side of the valley.

This is changing the rate of debris flow activity, and concentration in a specific area influences landslides susceptibility and magnitude–frequency. The change to the landslide hazard has an impact on the estimated risk levels. The increase in activity from a particular area is a common phenomenon, based on the long-term observations of national park staff and glacier guides (Marius Bron, personal communication, 2018), with this type of behaviour described colloquially as “switching on and off”, whereby a particular gully or slope will display enhanced rates of landslide activity for a period of time (sometimes on the order of years) before the levels of activity reduce. We incorporated the elevated debris flow activity into the risk analysis by deriving a specific magnitude–frequency relationship for debris flows from the Alpine Gardens and Mill's Creek catchment and applying this revised magnitude–frequency relationship to source areas within this catchment (see the Appendix for information).

The increased landsliding from the Alpine Gardens area was propagated through the risk equation. Figure 11 displays the LPR map that includes the current elevated rates of debris flow activity in the Fox Glacier valley. In the Alpine Gardens–Mill's Creek catchment, the increase in LPR ranges from 5 % to 1442 %, with a mean increase in LPR of 285.5 % ± 245 %.

Figure 11 Calculation of local personal risk for the Fox Glacier valley, including recent elevated levels of activity in the Alpine Gardens–Mill's Creek debris complex.

5 Discussion

5.1 Drivers of risk

We disaggregate our QRA to determine the dominant contributors to risk in each glacier valley. For both valleys, in our central-estimate scenario, aseismic landslides dominate the risk profile. Of these aseismic landslides, the major contributors to the risk are the moderately sized landslides (10⁴ to 10⁵ m³), which happen more frequently than the large or very large landslides (>10⁵ m³) but travel further and impact a larger area than the more frequently occurring small landslides (<10⁴ m³). Only when the tracks veer closer to the base of the slope does the risk from small landslides and rockfalls increase. We suggest that a similar pattern would be observed for seismic landslide volumes, given that the same volumes will impact the same area and that increases in the risk associated with the 10⁴ m³ landslide volume class in Franz Josef are associated with increases in both seismic and aseismic LPR.

Even when the annual frequency of a large ground-shaking event is increased to account for the probability of a >M_w 8 Alpine Fault event occurring in the next 50 years (Howarth et al., 2021), aseismic landslides account for more than half of the risk in the Fox Glacier valley. In the Franz Josef Glacier valley, increases in the number of aseismic landslides result in a large increase in risk of 365 %, suggesting that within the aleatory uncertainty (e.g. within the standard deviation) in the landslide magnitude–frequency relationships, aseismic landslides are a dominant contributor to the risk. Similar conclusions were reached by Robinson et al. (2016) in their analysis of coseismic landsliding from an Alpine Fault event, which suggests that for the central section of the Southern Alps, aseismic erosional processes are more important than seismically driven landslide erosion on annual timescales.

However, both accounting for the increased probability of an Alpine Fault earthquake occurring in the next 50 years and increasing the number of landslides that could occur during an earthquake event lead to increases in our risk estimates, both the individual risk per trip, discussed in detail here, and societal risk, as determined by f–N pairs, which represent the frequency (f) of an accident killing one or more people (N) in a single event. In Franz Josef, increasing earthquake annual frequency and the number of seismic landslides resulted in increases in risk of 330 % and 46 %, respectively. The increase in earthquake annual frequency results in a larger cumulative increase in risk of 330 % compared to 46 % for the increased number of seismic landslides – the difference may be the result of the order of the scenarios and compounding effect of variables. For example, the cumulative increase in risk associated with including a time-dependent earthquake scenario in Scenario 7 is 330 % and 56 % for Franz Josef and Fox, respectively, while for Scenario 9 (time-dependent-only scenario) it is 260 % and 54 %. We suggest that the differences between the example scenarios are due to changes in the number of seismic landslides generated, spatial probability of impact, and vulnerability. However, the relative differences between the scenarios do not change. In Fox, increasing the annual frequency of earthquakes and the increasing number of aseismic landslides result in risk increases of 56 % and 30 %, respectively. Increasing the number of seismic landslides results in risk increases of 10 %, lower risk increases than those observed in Franz Josef. In the Fox Glacier valley, the presence of large debris fans indicates debris flow activity (Gomez and Purdie, 2018; e.g. Cody et al., 2020); however, debris flow records in both valleys are limited, and therefore the aseismic debris flow risk may be underestimated. Our example from the Mill's Creek debris fan highlights that local increases in debris flow activity can significantly affect the risk, with local increases in risk of up to 1442 %. For both seismic and aseismic landslides, the impact of the number of landslides generated, which is the P_(L) term in the risk equation, emphasises the importance of the landslide inventory as an input into the risk calculation process. Therefore, more time and resources dedicated to the creation of a landslide inventory may reduce the uncertainty associated with the risk values (van Westen et al., 2008).

For seismic landslides, the landslide inventories of the 2016 M_w 7.8 Kaikōura, 1968 M_w 7.1 Inangahua, and 1929 M_w 7.8 Murchison earthquakes (Massey et al., 2018; Hancox et al., 2014, 2015) were used as proxies for the Franz Josef Glacier and Fox Glacier valleys given the lack of seismic landslide inventories for the west coast. All three inventories were dominated by shallow debris avalanches, with such failure types potentially being the dominant type of seismic landslide type (Keefer, 2002). The schist rock mass of both glacier valleys is fractured with persistent faulting (Cox and Barrell, 2007), and therefore we assume that shallow debris avalanches are the dominant failure type. While all three inventories occur in mountainous terrain similar to that of the Franz Josef Glacier and Fox Glacier valleys, climatic differences exist, with the impact of these climatic differences on the number and size of seismic landslides triggered unknown.

The biggest increase in risk values in Fox is associated with increases in the vulnerability and spatio-temporal probability of a visitor being in the path of a landslide (66 % and 80 % increases in risk, respectively), with these factors resulting in increases in risk in Franz Josef (219 % and 328 %, respectively). This emphasises the importance of risk management decisions to reduce exposure and lower vulnerability. Changes to the diameter of a person and debris height had a very limited impact on the estimated risk values, which affect the P_(T:L) term in the risk equation. Changes in the spatial extent of debris from the numerical simulations were not included within the sensitivity analysis but could be included using empirical or other probability-based (e.g. Flow-R – Horton et al., 2013) runout analysis and calculations of runout probability of exceedance (McDougall, 2017; e.g. Brideau et al., 2020). In this case study, we assume such variations are not particularly meaningful given the number of source areas the rockfall and landslide runout were simulated from (Fig. 2), the confined nature of the valleys, and the proximity of the access tracks and road to the base of the steep slopes from which the debris is sourced.

5.2 Time-variable risk

Our sensitivity analysis highlights the importance of accounting for time-variable risk with the inclusion of the increased frequency of an Alpine Fault earthquake resulting in cumulative increases in risk of 330 % and 56 % for the Franz Josef Glacier and Fox Glacier valleys, respectively. Alongside this, increases in aseismic landsliding result in cumulative increases in risk of 365 % and 30 % for the Franz Josef Glacier valley and Fox Glacier valley, respectively. Use of the upper estimate of the number of aseismic landslides that could occur may represent a future climate change scenario, reflecting the increased rates of landsliding (Gariano and Guzzetti, 2016) as glaciers retreat, slopes debuttress, and the environmental condition changes. Our sensitivity analysis suggests that climatically driven increases in landsliding will have a larger impact on landslide risk in the Franz Josef Glacier valley than in the Fox Glacier valley. We hypothesise that the differences in sensitivity analysis between the valleys may reflect the geomorphology of each valley. In Franz Josef, increases in the number of larger landslides may significantly increase the probability of such landslides reaching the element at risk.

However, the size and frequency of landslides may change in response to climate change (Huggel et al., 2012; Korup et al., 2012), with Liu et al. (2021) observing a shift in the frequency–area distribution with larger landslides occurring in their dataset of landslides in the high mountains of Asia. In our sensitivity analysis we do not test changes in the gradient of the magnitude–frequency distributions to reflect increases in the frequency of larger landslides occurring relative to the frequency of smaller landslides. Such shifts in the magnitude–frequency distribution will impact the risk results and associated uncertainty. The larger debris and alluvial fans in Fox may indicate higher rates of aseismic landsliding than those observed in Franz Josef and consequently may also explain that due to the already high rates of landsliding, increases in landslide rates (both seismic and aseismic) may have limited impact on risk, while the changes in vulnerability and exposure have a relatively bigger impact on the overall risk value.

Changes in landslide susceptibility should also be accounted for, as highlighted by Reichenbach et al. (2018), where the spatial pre-disposition to landsliding may change in response to environmental changes though the exact changes to landslide susceptibility are unknown. Given that landslide susceptibility is usually the starting point for risk analyses, time-variable landsliding and therefore susceptibility mean that the risk to people and infrastructure from landslides is also time-variable, especially after a major earthquake (e.g. Massey et al., 2014; Lin et al., 2006; Marc et al., 2015). Following the Canterbury earthquake sequence in New Zealand, a time-varying seismic hazard model was used as input to quantify the risk to life from rockfall in the Port Hills of Ōtautahi / Christchurch (see Massey et al., 2014). Consequently, the rockfall risk was shown to be time-variable with a rapid 50 % decrease in seismic rockfall risk in the 5 years post-earthquake event and a 14 % decrease in risk 5 to 10 years post-event. Massey et al. (2022b) show that aseismic rockfall risk is also elevated post-earthquake event with a similar 50 % decrease in rockfall rates 1 to 5 years post-earthquake event. Our example of increased debris flow activity on the Mill's Creek debris fan allows for spatial changes in landslide susceptibility frequency and magnitude to be easily incorporated into the risk model. Our analysis shows that the impact on LPR on the Mill's creek debris fan is significant. The ability to dynamically update the risk model to account for increased landslide activity in a specific area or catchment allows changes in environmental conditions and progressive failure, for example, increased number and size of landslides (Purdie et al., 2015; Fischer et al., 2012; Liu et al., 2021; e.g. Allen and Huggel, 2013) in a recently deglaciated area, to be assessed. The only required input is an estimate of approximate landslide size and frequency for a particular spatial area. Sensitivity analysis could be undertaken to understand if variations in the magnitude–frequency relationship had a significant impact on the resulting risk estimates. It is also important to note that our risk analysis does not include cascading hazards, such as landslide dam formation and associated dam break floods as well as catastrophic glacier multi-phase mass movements, which may be important in an Alpine Fault earthquake scenario (Robinson and Davies, 2013). Such cascading hazards could be incorporated into future risk analysis potentially using an event tree approach (e.g. Macciotta et al., 2016). Alongside changes in hazard behaviour, risk analysis should also account for dynamic changes in exposure and vulnerability. Voumard et al. (2013) developed a dynamic traffic simulator to simulate changing traffic speeds, with their results showcasing increased risk due to slow-moving traffic and traffic light placement compared with static speed and therefore increased exposure. The Stock et al. (2014) quantitative rockfall risk analysis of people within in the Yosemite Valley highlights that closure of buildings within their rockfall hazard line and reduced exposure resulted in decreased cumulative risk in the valley. Since 2019 and 2020 the main visitor tracks in the Fox Glacier and Franz Josef Glacier valleys, respectively, have been closed or partially closed due to geomorphic processes. Until access is restored in both valleys, the exact location of the tracks in each valley and the number of people walking the tracks are unknown. In the Fox Glacier valley, the expansion of the Mill's Creek fan from debris flow activity has damaged access on the true right side of the valley, while in the Franz Josef Glacier valley, the course of the Waiho River has restricted access on the true left side of the valley. As such visitor exposure and therefore risk to landslide hazard are reduced. Alongside this, the Covid-19 pandemic and associated closure of New Zealand's border to international tourists have resulted in a reduction in visitor numbers to both glacier valleys. This reduction in visitor numbers will impact our societal risk metric, by reducing exposure of one or more people to an event that might result in fatalities.

Figure 12 Quantitative estimates of individual risk per trip for Fox Glacier and Franz Josef Glacier visitors compared against popular tourist activities, modes of transportation used to access the glaciers, and individual risk per trip to national parks overseas (sourced from Taig, 2022a, b). The band range for each activity represents both the statistical uncertainty and uncertainty in the denominators of units of activity undertaken (see Taig, 2022a, b).

5.3 Risk communication and management

Our analysis quantitatively estimates the risk to life to visitors from landslides, with this information used by risk managers and decision makers to evaluate risk tolerability, determine appropriate risk mitigation measures, and communicate the risk to visitors and workers in each valley. Due to the uncertainty associated with risk analysis (Lee and Jones, 2014), we report our risk estimates as bands and not as single points (see Fig. 12). The risk bands represent our central estimate through to our upper estimate of the risk – we do not present a lower estimate of the risk as lower estimates are not currently used in decision-making regarding risk acceptability – in order to ensure that the highly uncertain risk levels are not underestimated. However, we note and agree with Strouth and Mcdougall (2021), who state that risk assessment conservatism should be avoided, with central estimates used for risk evaluation and uncertainties in the risk analysis presented transparently. The risk bands can be presented against risk comparator data to inform risk evaluation and risk tolerability processes in conjunction with an evaluation of how visitors and decision makers perceive risk (see Taig et al., 2012; Taig, 2022a, b). In Fig. 10, the individual risk per trip for visitors to both the Franz Josef Glacier valley and the Fox Glacier valley (though it is important to note that the risk numbers do not include any mitigation measures and therefore do not represent residual risk) is plotted against other activities that a visitor may undertake. These activities include popular tourist activities in New Zealand, modes of transport to and from the glacier valleys, and risk per trip in other national park settings globally. More information on these datasets can be found in Taig (2022a, b). Not only can this be used to inform the risk evaluation process, but it can also help with risk communication to visitors. The range in risk values can be presented as a graphic to illustrate the risk and to avoid confusion with small numbers or scientific notation, as well as to help visitors, whose main language may not be English, understand the uncertainty in risk results (Taig, 2022a, b).

The disaggregation of the QRA allows a greater understanding of both the contributors to landslide risk and their associated uncertainty. Such an approach presents a useful tool to inform and communicate to risk managers where appropriate management and mitigation strategies may be most effective. Reductions in vulnerability and exposure can be important risk mitigation measures (e.g. Schneiderbauer et al., 2017), as highlighted by our risk sensitivity scenario analysis. The LPR maps can be used to inform track placement and realignment, reducing an individual's time spent and exposure in high-hazard zones beneath steep slopes. The LPR maps may also be used for identifying suitable, “less risky” stopping points in the valley, when tracks are partially closed due to rainfall, high streamflows, or other events. Where it is not possible to relocate the tracks, other mitigation measures such as track closures during heavy rainfall may reduce the risk from rainfall-induced landslides within the aseismic landslides class. As aseismic landslides dominate the risk profiles, reduction in exposure to rainfall-induced landslides may result in a significant reduction in visitor risk per trip. In Fig. 12, we include a theoretical reduction in aseismic landside risk of 75 %, assuming that 75 % of aseismic landslides are triggered under heavy-rainfall conditions. In this example, if the track is closed under heavy-rainfall conditions, substantial reductions in the visitor risk per trip occur (see Massey et al., 2022a). However, due to the limitations of our landslide inventory, we are unable to link landslide occurrence to rainfall events in each valley and therefore cannot provide a quantitative basis for the risk reduction associated with track closures. To provide a robust basis for using track closures as a risk reduction method, rockfall and landslide events should be documented and recorded within each valley, along with meteorological observations. This should also include any information on the occurrence of debris flows, particularly within Fox, as this landslide type is difficult to determine within our landslide inventory analysis and is therefore underrepresented in our magnitude–frequency analysis even though debris fans are higher-risk environments (see Fig. 6). Such information could be crowdsourced, similarly to in the Yosemite Valley, for example, where visitors are able to report rockfall occurrences to park staff (https://www.nps.gov/yose/learn/nature/rockfall.htm, last access: 6 July 2022). For both the Franz Josef valley and the Fox Glacier valley, the glacier guides continue to record rockfall and landslide activity. A rockfall/landslide register can also be used to inform dynamic risk analysis by recording areas of locally high activity. Our methodology presents a base risk model that can be easily updated and amended to incorporate future information such as revised track locations, visitor numbers, and changes in landslide activity.

6 Conclusion

We have presented a quantitative risk analysis (QRA) case study from the Franz Josef Glacier and Fox Glacier valleys, on the west coast of the South Island, New Zealand. We deconstructed the QRA to reveal the relative contributions of aseismic versus seismic landsliding and landslide volume classes to risk. Our results reveal that for both valleys in our central-estimate scenario, aseismic landslides contribute more to the overall risk than seismic landslides. However, our sensitivity analysis of nine risk scenarios, to explore the uncertainties in our inputs to the model, suggests that the contribution of seismic or aseismic landslide risk is dependent on time-variable input assumptions. The increasing probability of a large Alpine Fault earthquake occurring results in increased seismic landslide risk, both individual and societal. Increases in the number of aseismic landslides, within the standard deviation of the valley-specific magnitude–frequency relationships, also increase the landslide risk, particularly in Franz Josef. This increase in aseismic landsliding may reflect climatically induced changes in landslide rates in these actively deglaciating valleys and suggests that the risk of landsliding will change under different climate change scenarios. Additionally, the spatial location and susceptibility of landsliding may also change in response to environmental changes. We presented an example to show-case how local changes in the rates of landsliding can be explored and incorporated in the analysis. We present our risk results as bands, not points, that display the uncertainty in our risk results. We suggest that QRA not only is a valuable tool for evaluating the risk to an individual but can be used to better understand what drives landslide risk and as such what risk management decisions will be most effective and appropriate in significantly reducing risk. In order to do this, QRA must be able to be deconstructed as well as be dynamic to account for changing hazards and exposure with time.

Appendix A: Aseismic landslide inventories

The data sources for the aseismic landslide include (1) a rockfall register compiled from observations made by staff of the Department of Conservation (DOC), Franz Josef Glacier Guides Ltd, and Fox Glacier Guiding; (2) a landslide inventory derived from historical aerial imagery analysis of both valleys; and (3) a large landslide inventory of historical landslides observed in the wider Southern Alps. The record of observed rockfall activity in the rockfall register contains data collected since 2008 for Franz Josef and 2009 for Fox. The rockfall registers record the date, approximate size, and source location, if identifiable, of rockfalls. Alongside this, local knowledge of long-term guides Craig Buckland and Jon Tyler (Franz Josef Glacier Guides Ltd; personal communication, 2018) and Marius Bron (Fox Glacier Guiding; personal communication, 2018) informed the relative changing rates and sources of landslide activity within each valley. We identified and mapped landslides from a series of historical aerial photographs for each valley (1948, 1965, 1981, 1985, 1987, 2011, and 2017 for the Franz Josef Glacier valley, and 1953, 1981, and 2017 for the Fox Glacier valley), and these landslides were subsequently verified in the field. Additionally, we identified large ($>\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ m³) relict landslides in each valley, with an unknown temporal occurrence. To assess the potential for large landslides to occur, we used a landslide dataset recorded in the wider Southern Alps region, where there is evidence of large landslides occurring under aseismic conditions, such as the 11×10⁶ m³ Young River landslide in 2007 (Massey et al., 2013). We also used data from the following studies, which detail the occurrence of debris avalanches since 1978: McSaveney (2002), Hancox et al. (2005), Cox and Allen (2009), Allen et al. (2011), Allen and Huggel (2013), Massey et al. (2013), and Cox et al. (2015). We assumed that landslides $<\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ m³ are unlikely to have been noticed or mapped unless they impacted people or property in the wider Southern Alps. We also assume that landslides in both valleys where the glacier guides and DOC operate would be well documented as people are present almost on a daily basis.

Appendix B: Aseismic landslide susceptibility models

We used the best sub-set regressions to explore which group of variables could statistically best explain landslide occurrence. From these variable groupings, we undertook backward stepwise regression modelling to determine which group of variables was the most statistically significant. Using the variables of slope angle, local slope relief (LSR), material type, and vegetation, we estimated the aseismic landslide probability using the following logistic regression equation (with output coefficients in Table B1) for the Fox study area:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& \mathrm{non}\text{-}\mathrm{EQ}\mathrm{LS}\mathrm{probability}\\ & =\mathrm{1}/\left(\mathrm{1}+e^{\left(-\left(\mathrm{intercept}+\mathrm{slope}\mathrm{angle}+\mathrm{LSR}+\mathrm{material}\mathrm{type}+\mathrm{veg}\right)\right)}\right).\end{array}\end{array}$$

For the Franz Josef study area, we found vegetation to be a statistically insignificant variable when used to explain landslide occurrence, and as such it was not included in the model. For the aseismic landslide probability in the Franz Josef study area, we used the following logistic regression equation (with output coefficients in Table B2):

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}& \mathrm{non}\text{-}\mathrm{EQ}\mathrm{LS}\mathrm{probability}\\ & =\mathrm{1}/\left(\mathrm{1}+e^{\left(-\left(\mathrm{intercept}+\mathrm{slope}\mathrm{angle}+\mathrm{LSR}+\mathrm{material}\mathrm{type}\right)\right)}\right).\end{array}\end{array}$$

Table B1 Summary table of coefficient estimates for the variables used in the Fox study area logistic regression equation.

^a Colluvium is set as the reference material, which means that alluvium is less likely to fail with a negative estimate and rock is more likely to fail with a positive estimate. ^b Slopes that are not vegetated are set as the reference vegetation variable, which means that areas that are vegetated are less likely to fail with a negative estimate. n/a: not applicable.

Table B2 Summary table of coefficient estimates for the variables used in the Franz Josef study area logistic regression equation.

^* Alluvium is set as the reference material, which means that colluvium and rock are more likely to fail with a positive estimate. n/a: not applicable.

Appendix C: Landslide runout analysis

C1 Rockfall

We modelled landslides with volumes ≤ 1000 m³ as rockfalls using RAMMS rockfall software (RAMMS, 2015) for all material types. The software simulates the rigid-body motion of falling rocks and predicts rock trajectories in general three-dimensional terrain. Rock trajectories are governed by the interaction between the rock, its associated shape, and the nature of the ground (e.g. a soft substrate such as sand will dampen the rock energy in contrast to a hard substrate such as rock). Generalised rock shapes are simulated, and rock block orientation and rotational speed are included in the rock–ground interaction. We determined the simulation parameters for forecasting by back analysing recorded rockfalls within the study areas, where the source area, boulder shape, and rockfall trails were recorded or could be accurately inferred. The RAMMS rockfall forecast parameters adopted from back analysis are shown in Table C1, along with descriptions of the parameters and the data sources used to derive them. The results from the simulations comprise the following: kinetic energy, runout distance, jump heights, and the number of simulated trajectories passing through a given grid cell.

C2 Debris avalanche and flows

To identify the areas impacted by landslides and the associated height of debris and number of boulders, we conducted a suite of runout simulations for the different volume bins. We modelled landslides with source volumes > 1000 m³ as debris avalanches (if sourced from rock) or as debris flows (if sourced from colluvium or moraine) using RAMMS debris flow software (RAMMS, 2011). RAMMS is based on the Voellmy friction law, where the frictional resistance consists of a dry-Coulomb-type friction (coefficient μ), which scales with normal stress, and a viscous turbulent friction (coefficient ξ), which scales with landslide volume. These coefficients are calibrated from the back analysis of case studies. For this assessment, we used the back analysis of 67 debris avalanches (ranging in volume from 300 to 100×10⁶ m³) published in the literature (Schneider et al., 2011; Allen et al., 2009). For debris flows, we used 22 back analysis case studies ranging in volume from 1000 to 2×10⁵ m³ (Loup et al., 2012; Cesca and Agostino, 2008; Deubelbeiss et al., 2011; Hussin, 2011; Scheuner et al., 2011). We fitted a power law to the data (Figs. C1 and C2) to calculate the coefficients for the numerical simulations. For debris flows, the ξ parameters did not vary with source volume, and so we adopted a central estimate of 350 in the numerical simulations.

In areas where the source area could potentially fail as either a debris avalanche or debris flow (for example, potential failure from the top of the larger creeping Yellow Creek landslide in the Fox study area), we simulated both and calculated the maximum debris flow height from the two outputs.

C3 Sensitivity analysis of debris avalanche and debris flow output

We assessed the sensitivity of the simulated maximum debris heights to varying RAMMS input parameters. Different input μ and ξ parameters within RAMMS result in a change in both the extent of the debris runout (and therefore area inundated by debris) and the height of the debris and therefore the number of boulders passing through a given location on the ground (grid cell). For the debris avalanche simulations, we calculated the standard error of the modelled fit of the data (Fig. C3). We both added and subtracted the standard error from the power-law relationship to obtain the mean ± 1 SE (standard error) values of both μ and ξ parameters (Fig. C3). For the debris flow simulations, we used the same procedure as outlined for debris avalanches to calculate the μ parameter. As no relationship existed for the ξ debris flow parameters (Fig. C2), we used the standard deviation (σ) of the mean parameters to calculate both mean + 1σ and mean − 1σ (Fig. C4). We choose one representative source area for both debris flow and debris avalanche deposits to simulate both mean + 1 SE (or standard deviation) and mean − 1 SE (or standard deviation) runout. For each source area, the following volume classes were simulated; 10⁴, 5×10⁴, 10⁵, and 10⁶ m³. We varied the simulated volume size to assess if the range in maximum flow height increased or decreased for larger volumes. From the simulation results, we calculated the difference in debris heights (per grid cell) between the standard parameter simulation results and the results from each respective mean ± 1 SE simulation. We calculated the mean and standard deviation of the difference (in debris height per grid cell) for each simulation result (i.e. we calculated the difference in the difference). We summed the values to calculate the mean + 1σ (upper-bound) value of the difference between the simulations.

The results from our sensitivity assessment indicate that the absolute difference in debris heights increases with volume size, whereby larger landslides can display several metres of difference in flow height for any given grid cell. Proportionally, the differences in maximum flow height for debris avalanches were on average 60 % ± 22 % higher than those modelled using the preferred forecast parameters. For debris flows, the difference in maximum flow heights were on average 60 % ± 28 % higher than those simulations adopting the preferred forecast parameters. We applied this 60 % factor of difference to all simulation results to derive upper estimates of debris heights.

Table C1 RAMMS rockfall model parameters used for forecasting rockfalls.

n/a: not applicable.

Figure C1 Debris avalanches: range of parameters used to back analyse the runout of debris avalanches published in the literature (n=67), using the RAMMS software.

Figure C2 Debris flow: range of parameters used to back analyse the runout of debris flows published in the literature (n=22), using the RAMMS software.

Figure C3 Range of parameters used to back analyse the runout of debris avalanches published in the literature. The fitted purple line represents the mean − 1σ; the red line represents the modelled fit of the data; the green line represents the mean + 1σ.

Figure C4 Range of parameters used to back analyse the runout of debris flows published in the literature. The fitted purple line represents the mean − 1σ; the red line represents the rounded parameters used within the simulation; the green line represents the mean + 1σ.

Appendix D: Calculation of LPR

For seismic landslides, the probability of death (P_D) is calculated for each earthquake Band 1 to 4, for each individual source area of a given volume class and the debris that the given source area generates. The calculations were done for each grid cell within the study areas. Firstly, we calculated the probability of death for each landslide source and its related debris (P_D(Source)) as follows:

$$\begin{array}{}\text{(D1)}& \begin{array}{rl}& P_{\mathrm{D}(\mathrm{EQ}\mathrm{band}x;\mathrm{vol}\mathrm{class}y;\mathrm{source}z)}\\ & =P_{A\mathrm{or}B\left(\mathrm{dependent}\mathrm{on}\mathrm{LS}\mathrm{vol}\mathrm{class}\right)}\times P_{\mathrm{2}(\mathrm{S}:\mathrm{H})}\times V.\end{array}\end{array}$$

P_A or B is the probability the given source within each landslide volume class will generate the given volume of debris. P₂ is the probability of being in the path of N boulders within the debris generated by a given source if it occurs. V is the vulnerability, defined as the probability of a person being killed if present and in the path of one or more boulders. We then calculated the P_D(source) for each source area of a given landslide volume class (y) and its related debris and combined this for each landslide volume class to estimate the probability of death from all landslides of the same volume class that might contribute to the given grid cell (P_{D(vol class)}):

$$\begin{array}{}\text{(D2)}& \begin{array}{rl}& P_{\mathrm{D}(\mathrm{EQ}\mathrm{band}x;\mathrm{vol}\mathrm{class}y)}\\ & =\mathrm{1}-\left(\mathrm{1}-P_{\mathrm{D}\left(\mathrm{vol}\mathrm{class}{y}_{\mathrm{1}};\mathrm{LS}\mathrm{source}{z}_{\mathrm{1}}\right)}\right)\\ & \times \left(\mathrm{1}-P_{\mathrm{D}\left(\mathrm{vol}\mathrm{class}{y}_{\mathrm{1}};\mathrm{LS}\mathrm{source}{z}_{\mathrm{2}}\right)}\right)\\ & \times \left(\mathrm{1}-P_{\mathrm{D}\left(\mathrm{vol}\mathrm{class}{y}_{\mathrm{1}};\mathrm{LS}\mathrm{source}{z}_N\right)}\right).\end{array}\end{array}$$

We then calculated the probability of death from all landslides of a given volume class generated by each earthquake band. P_{D(EQ band x; vol class y)} for a given earthquake band (x) and volume class (y) was multiplied by the number of landslides of a given volume class (N_LS)) generated by the representative earthquake PGA of the given band. If the number of landslides (N_LS) triggered in the band was ≥1, then instead of multiplying P_{D(EQ band x; vol class y)} by the number of landslides (N_LS), the following formula was used:

$$\begin{array}{}\text{(D3)}& \begin{array}{rl}& P_{\mathrm{D}(\mathrm{EQ}\mathrm{band}x;\mathrm{all}\mathrm{vol}\mathrm{class}y)}\\ & =\mathrm{1}-{\left(\mathrm{1}-P_{\mathrm{D}\left(\mathrm{EQ}\mathrm{band}x;\mathrm{vol}\mathrm{class}y\right)}\right)}^{N_{\mathrm{LS}}}.\end{array}\end{array}$$

We combined the contribution to each grid cell from each landslide volume class per band to calculate the probability of death from all landslides triggered by the representative PGA in the given band (P_{D(EQ band x)}):

$$\begin{array}{}\text{(D4)}& \begin{array}{rl}& P_{\mathrm{D}\left(\mathrm{EQ}\mathrm{band}\mathrm{x}\right)}=\mathrm{1}-\left(\mathrm{1}-P_{\mathrm{D};\mathrm{all}\mathrm{vol}\mathrm{class}\mathrm{1}k}\right)\\ & \times \left(\mathrm{1}-P_{\mathrm{D};\mathrm{all}\mathrm{vol}\mathrm{class}\mathrm{10}k}\right)\times \left(\mathrm{1}-P_{\mathrm{D};\mathrm{all}\mathrm{vol}\mathrm{class}N}\right).\end{array}\end{array}$$

We calculated the local personal risk from all landslides that occur within the given band by multiplying P_{D(EQ band x)} by the annual frequency of the representative earthquake PGA in that band.

For aseismic landslides, we calculated the probability of death (P_{D(vol class y)}) in the same way as for earthquakes except ignored the need to calculate for each earthquake band. We then calculated the LPR for each volume class by multiplying P_{D(vol class y)} by the annual frequency of the given volume class (y) of landslide occurring.

Appendix E: Mill's creek catchment magnitude–frequency relationship

We used several datasets to derive the magnitude–frequency relationship for the Mill's Creek catchment, including (1) a change detection model from differencing of a March 2017 digital surface model (DSM) and June 2018 digital elevation model (DEM); (2) national park staff observations of the frequency of debris flow events and a rough estimate of their associated volume; and (3) NIWA weather observation data from the village of Franz Josef, which represents the closest meteorological observation point. Our change detection model revealed that between March 2017 and June 2018, approximately 6.5×10⁶ m³ was eroded and 3×10⁶ m³ deposited within the Alpine Gardens and Mill's Creek catchment. During this same time period, the valley was closed 34 times due to heavy rain and flooding. For the larger storm events, including ex-tropical Cyclone Fehi in February 2018, national park staff observed debris flow activity that resulted in damage to the road. The staff estimated that for the Cyclone Fehi event, approximately 2×10⁶ m³ had been deposited on the Mill's Creek debris fan (Tony Hart, personal communication, 2018). Using this information, including both the rough national park staff volume estimates and the frequency of heavy-rain events likely to trigger debris flows, we divided the approximately 6.5×10⁶ m³ into different debris flow events based on magnitude–frequency principles. We normalised the data over the 1.38-year time record and the spatial area of the Alpine Gardens and Mills Creek catchment to derive a magnitude–frequency–power-law relationship, with elevated rates of landslide activity compared to the magnitude–frequency relationship for the Fox Glacier valley overall.