Introduction

Landslides are the second most frequent tsunami source worldwide (Harbitz et al. 2014a; Tappin 2010; Yavari-Ramshe and Ataie-Ashtiani 2016). Their importance is undisputed, most recently demonstrated by the 2018 Anak Krakatoa volcano flank collapse (Grilli et al. 2019) causing several hundred fatalities. Additional examples of fatal landslide tsunamis comprise among others the 1959 Vajont event (Crosta et al. 2016), the 1998 Papua Guinea tsunami (Tappin et al. 2008) and the 1792 Mount Unzen tsunami (Sassa et al. 2016). Moreover, several subaerial landslide-generated waves have induced up to 150-m run-up heights during the last two decades (George et al. 2017; Gylfadóttir et al. 2017; Paris et al. 2019; Sepúlveda and Serey 2009; Tinti et al. 2005). Further examples can be found in the review of (Harbitz et al. 2014a).

Landslide tsunamis have the potential to produce larger, but often more local (Okal and Synolakis 2004) tsunamis than earthquakes, and they are also more complex and display a more diverse nature (Løvholt et al. 2015). Therefore, prognostic analysis of landslide tsunamis is far less developed than for earthquakes, and hence, we have a more limited understanding of landslide tsunami hazards (Harbitz et al. 2014a; Geist and Parsons 2014). A first reason for this is due to lack of knowledge related to the landslide occurrence frequency, with just a few limited records of comprehensive landslide statistics covering landslide volumes across several orders of magnitude (Blikra et al. 2005; Geist and ten Brink 2019; Lane et al. 2016; Urgeles and Camerlenghi 2013). Hence, the statistics describing the landslide recurrence is more poorly constrained than for earthquakes and sometimes even non-existing. A second reason is that in real situations, we lack well-established procedures for estimating the uncertainty range the landslide dynamics would take. This concerns uncertainties related to estimating the material properties in situ, corresponding uncertainties related to geology and failure planes, as well as the complex material behaviour and transformation after the landslide have been released. So while we do have significant knowledge about how tsunamis are generated once the motion of the landslide is known, see, e.g. (Løvholt et al. 2015), the landslide dynamics and probability of release is not known a priori for potential hazardous events. This precise information about the landslide motion and release probability is difficult to obtain due to the enormous geographical scale that often needs to be taken into account, combined with a lack of relevant data. All these factors play major roles in quantifying the corresponding tsunami generation uncertainty and, hence, on how they propagate into corresponding tsunami uncertainties.

Consequently, past prognostic analyses of landslide tsunami hazard have historically resorted to deterministic analysis, e.g. (Harbitz et al. 2014b; Li et al. 2019; Schambach et al. 2019). Similar studies addressing the sensitivity to landslide parameters are more limited (Hilbe and Anselmetti 2015; Kim et al. 2019; Ren et al. 2019; Salmanidou et al. 2017). To this end, just a handful of landslide probabilistic tsunami hazard analysis (PTHA or LPTHA) studies also quantifying the temporal probability of tsunamis exists in the scientific literature, e.g. Grezio et al. 2012, Lane et al. 2016 and Grezio et al. 2020. However, none of the landslide PTHA studies to date appropriately addresses tsunami uncertainty due to the unknown variability in landslide dynamics. There is hence a pressing need to better understand and prepare for landslide tsunami hazard. This encompasses both the understanding of the landslide tsunami uncertainty, understanding their frequency of occurrence and setting up a suitable probabilistic framework to carry out the analysis. In this paper, we attempt to shed light on these processes and the computational framework by:

  1. 1.

    Providing a review of recent studies that quantify landslide dynamics uncertainty and its effect on tsunami generation, presented in the “Landslide dynamics uncertainty constrained by data: results from previous studies” section.

  2. 2.

    Presenting a simple event tree based on a probabilistic framework for PTHA, incorporating different uncertain parameters related to landslide dynamics. Acknowledging that the framework is simplified, we argue that it still poses a sufficiently sophisticated structure for better addressing landslide dynamic uncertainty in PTHA than other models available today. This framework is presented in the “Framework for probabilistic tsunami hazard analysis: Lyngen case study” section.

  3. 3.

    Presenting case studies exemplifying how a probabilistic framework can be utilised in practical applications. The example case is the Lyngen fjord in Northern Norway, presented in the “Probabilistic example studies for Lyngen fjord” section.

  4. 4.

    In light of the findings in the technical sections sections, finally discussing in the “Concluding remarks” section how the landslide hazard models should be improved in the future, for instance to more rigorously take into account uncertainties.

Landslide dynamics uncertainty constrained by data: results from previous studies

Landslide tsunami genesis: a brief summary of past reviews

Modelling of landslide tsunami generation must take into account the time dependency of the landslide motion and can involve transitions from supercritical flow (the landslide outrunning the tsunami) to subcritical flow (the tsunami outrunning the landslide). In this respect, the ratio between the landslide velocity and the tsunami wave speed, defined as the Froude number Fr, measures the degree of criticality. This is in contrast to earthquakes, which most often rupture much faster than the tsunami propagates, meaning that the time dependency is relatively less important for the tsunami hazard. The landslide dynamics hence adds additional complexity to the tsunami hazard analysis that can be simplified for earthquakes. To this end, a range of factors controls the tsunami genesis, as discussed extensively in the literature, e.g. (Yavari-Ramshe and Ataie-Ashtiani 2016; Løvholt et al. 2015; Fritz et al. 2003; Grilli and Watts 2005; Harbitz et al. 2006; Haugen et al. 2005; Heller et al. 2008; Løvholt et al. 2005; Ruffini et al. 2019; Ward 2001; Watts 2000; Zengaffinen et al. 2020), and we only provide a simplified review here, presenting a few essential findings:

  • Different types of landslides can have extremely different tsunami-genic properties. The most pronounced difference is between subaerial landslides and submarine landslides.

  • The tsunami genesis is primarily controlled by the landslide volume and for impact tsunamis, frontal area, while the landslide dynamics (velocity and Froude number, acceleration) also plays an important role.

  • For subaerial landslides, the landslide tsunami genesis is studied extensively through a range of laboratory studies. Empirical analysis reveals that the generation scales non-linearly with the products Froude number and the frontal slide area cite Fritz et al. (2003), while later studies have factored in the effect of the landslide mass and slope angle (Heller et al. 2008; Bullard et al. 2019), the landslide length (Enet and Grilli 2007; Panizzo et al. 2005) and 3D directivity effects (Evers et al. 2019; Heller and Spinneken 2015; Mohammed and Fritz 2012).

  • For submarine landslides, the scaling can still depend on the Froude number, but is more strongly dependent on the acceleration (Haugen et al. 2005; Løvholt et al. 2005; Watts 2000; Hammack 1973) and has a different Froude scaling relationship (Løvholt et al. 2015; Ward 2001) compared with subaerial landslides. For submarine landslides, the tsunami genesis is then most efficient when Fr = 1, when the landslide and tsunami travel at the same speed.

While the landslide dynamics is uncertain, two main avenues can be considered for retrieving more accurate information about the tsunami genesis: The first one is using data from landslides and landslide tsunamis to quantify the landslide tsunami genesis. The second one is to combine geotechnical data with models for the landslide dynamics. In the remaining part of this section, we will highlight how the first approach, using landslide data and associated tsunami observations, can better shed light on relevant landslide dynamics and tsunami uncertainties.

Examples employing landslide run-out data only

Here, we review briefly two case studies conducted by other research groups. The first example (Fischer et al. 2015) concerns a study attempting to invert a snow avalanche deposit into bounded probability distributions for the associated parameters for the slide dynamics. This study used a frictional-collisional avalanche model, computing uncertainty bounds for the landslide friction parameters by optimising the agreement with landslide outputs. While snow avalanches are different from our present application, the frictional-collisional model is very similar to submarine landslide models and hence relevant. The first iteration of their study (Fischer et al. 2015) showed that landslide uncertainties were excessive, and for some parameters, it was even difficult to properly constrain the uncertainty. A more recent update by some of the same authors (Kofler et al. 2017), using the Markov Chain Monte Carlo (MCMC) method, yielded smaller uncertainties for the same avalanche parameters. In conclusion, the authors found that it was possible to obtain estimates of the most probable snow avalanche input parameters, but that the related uncertainty bounds were rather large.

The second example concerns ensemble simulations mimicking a large submarine landslide at the Rockall Bank offshore Ireland (Salmanidou et al. 2017). The procedure here was very similar to the one above, the frictional-collisional landslide model Volcflow Kelfoun et al. 2010 was used for the landslide simulation realisations, and MCMC was used to optimise the landslide parameters. A range of different maximum velocities was obtained, spanning the wide range of 30–130 m/s. One essential parameter for determining the velocity, the turbulence coefficient, exhibited an almost unitary probability distribution. Hence, landslide dynamics uncertainties were rather large. Other parameters, such as the basal friction and the yield strength, were more confined. When fed into a tsunami model, these uncertainties resulted in tsunami surface elevations ranging from 5–25 m for 100 random realisations. It is stressed that this uncertainty range is related to the tsunami generated by a single landslide volume at a given coastal point.

These two example cases indicate that the presence of landslide data alone can be used to estimate landslide uncertainties, but that large uncertainties should be expected, even when the volume is given. In the next couple of subsections, we will demonstrate how tsunami data can drastically reduce these uncertainties.

Subaerial landslide tsunami case: 2014 Lake Askja

A 20-Mm3 landslide impacted Lake Askja in Iceland on 21 July 2014 (Gylfadóttir et al. 2017). The landslide caused a tsunami with run-up heights ranging from 10–80 m. Because the event impacted a lake of limited size, both the landslide volume, run-out distance and run-up heights, were thoroughly documented. In particular, the tsunami run-up height was relatively densely sampled along the entire periphery of the lake. Because of the thorough documentation of the event (arguably the best-documented historical landslide tsunami to date), combined with the limited size of the lake, the Lake Askja event is viewed as ideal for studying and understanding landslide tsunami genesis.

The Lake Askja landslide tsunami was modelled in detail by Gylfadóttir et al. (2017). They used a simple landslide block slide model assuming a linear Coulomb term for the bed resistance and viscous drag and added mass resistance forces due to the ambient water. This landslide block was in turn coupled to both a Boussinesq-type (Kim et al. 2017) and a shallow water-type tsunami propagation and inundation model Geoclaw (Berger et al. 2011) implemented in the Finite Volume solver Clawpack (LeVeque 2002). Only the simulations including dispersion were able to closely match landslide run-out and tsunami run-up observations.

Because there exists no measurement or eyewitness observation of the landslide dynamics, it was necessary to test different realisations of the block model kinematics. Several hundred simulations were conducted, and the results were compared with the tsunami run-up height distributions. To quantify the degree of fit, an objective function summing the squared differences between simulated and observed run-up heights was used. Using an optimization procedure, it was found that the most likely landslide impact velocity was 30.9 m/s, and correspondingly, the most likely landslide thickness during impact was 35.5 m.

In Fig. 1, we show histograms as functions of the impact velocity of the slide and the slide thickness displaying the uncertainty of the impact velocity derived from Gylfadóttir et al. (2017). We here introduce the objective function of fix referencing as objAsk, noting that small values of objAsk implies a good match with observations. In Fig. 1 only results for objAsk < 1.5 × min(objAsk) are shown. To this end, min(objAsk) is the realisation that gives the closest agreement with observations. It is argued that all these cases with small values of objAsk represent a good match with the field observations (see also Gylfadóttir et al., 2017). Evident from Fig. 1, we see that there is a significant spread in the landslide impact speeds with 30–35 m/s, while the distribution for the landslide thickness is more confined between 34 and 37 m. Nevertheless, we see that the availability of tsunami data leads to a much smaller uncertainty in landslide speed, compared with the landslide dynamics for the Rockall Bank discussed above, where no tsunami data were available. This indicates that the tsunami, as a footprint of the landslide dynamics, can be used to provide much more accurate information on landslide kinematics than only the landslide runout observation. It is stressed again that these uncertainties are the only representative for hindcasting the dynamics of an event that has already happened. On the other hand, such information of past events should be considered useful for prognostic calculations in the absence of other data.

Fig. 1
figure 1

Results from an attempt to optimise landslide dynamics by hindcasting the well-documented Lake Askja landslide and tsunami (Gylfadóttir et al. 2017). The figure shows histograms for all the lowest values of the optimization parameter (objAsk< 1.5× min (objAsk)). Upper panel: distribution as a function of the landslide velocity at the moment of landslide impact with the water. Lower panel: distribution as a function of the landslide thickness

Full size image

Submarine landslide tsunami case: the Storegga Slide 8150 BP

The Storegga Slide (Bryn et al. 2005; Haflidason et al. 2004; Kvalstad et al. 2005) occurred 8150 BP offshore western Norway, and with a volume of about 3000 km3, it is one of the largest documented submarine landslides worldwide. Storegga resulted in an oceanwide tsunami, with palaeotsunami observations in Norway, Scotland, Faroe Islands and Denmark (Bondevik et al. 2005; Fruergaard et al. 2015; Løvholt et al. 2017; Smith et al. 2004). As for the Lake Askja rock slide, the relatively rich availability of both landslide and tsunami data for Storegga facilitates extensive testing of the landslide tsunami model. However, the much larger Storegga Slide involved a more complex landslide process, involving for instance retrogressive failure and material remoulding. These different mechanisms necessitated other kinds of modelling techniques. It was therefore necessary to use a more sophisticated landslide model, BingClaw (Kim et al. 2019; Løvholt et al. 2017), for modelling the Storegga Slide. BingClaw takes into account landslide remoulding to simulate the landslide dynamics. Kim et al. (2019) coupled BingClaw to the Geoclaw tsunami model, using the Storegga Slide to better understand the landslide dynamics sensitivity. In contrast to the Lake Askja study, the main focus of Kim et al. (2019) was to provide rough uncertainty bounds for the material parameters governing the slide motion, rather than just for the landslide kinematics as was done for Lake Askja. We briefly describe their main findings below.

Kim et al. (2019) studied the sensitivity to how the yield strength of the viscoplastic landslide material influenced by the landslide run-out distance and the tsunami run-up height. For the limited number of 64 case studies, they compared the simulated landslide runout with the observed run-out distance of the Storegga Slide, and with observed paleotsunami run-up heights. Yield strength remoulding between an initial strength τy,0 and a remoulded strength τy,∞ was assumed, assuming an exponential rate e−Γ γ where Γ is the rate factor and γ is the modelled material shear strain.

A roughly reworked summary of the findings of Kim et al. (2019) is presented in Table 1. The first row in Table 1 shows our interpreted uncertainty ranges based on the findings of Kim et al. (2019). It is stressed that this study did not undertake a parameter optimization study as described for Lake Askja in the “Subaerial landslide tsunami case: 2014 Lake Askja” section. Hence, the parameter ranges presented in Table 1 are subjective. The second row indicates which phases that are most sensitive to the given parameter according to Kim et al. (2019). Table 1 indicates that the landslide run-out distance is controlled by the remoulded yield strength of the material. On the other hand, Table 1 also indicates that the tsunami run-up height is controlled by the initial yield strength. This is supported by the previous studies (reviewed above) indicating that the initial acceleration phase of the landslide motion is controlling the tsunami genesis, and the initial yield strength strongly influences this acceleration phase. The rate at which the remoulding takes place has an influence on both the landslide run-out distance and the tsunami genesis.

Table 1 Rough material parameter uncertainty ranges (Kim et al. 2019) for BingClaw parameters applied for the Storegga Slide and tsunami
Full size table

While there are altogether relatively narrow uncertainty ranges for the initial and remoulded yield strengths, Kim et al. (2019) found that different parameter combinations could provide a relatively similar degree of match with paleotsunami run-up observations. The fact that different landslide properties influence different phases of the coupled landslide-tsunami process adds to the complexity. Hence, the multiplicity of landslide parameters controlling different phases of the landslide motion may lead to ambiguous interpretations of landslide parameter uncertainties when we use them in a forecasting situation. This uncertain must be interpreted as an epistemic uncertainty in prognostic modelling. On the other hand, the present parameter ranges are again more confined than the cases without tsunami observations reviewed above. As for Lake Askja, this demonstrates how tsunami information can be used to provide more accurate information on the landslide dynamics of past events.

Framework for probabilistic tsunami hazard analysis: Lyngen case study

Past, well-documented landslides and tsunamis provide relevant examples for constraining tsunami-genic landslide parameters. However, in most forecasting situations, we cannot expect that input landslide parameters and related probabilities can be completely constrained by field data. Hence, we need to resort to expert judgement in addition to empirical data, and the landslide parameters governing the tsunami genesis can be expected to exhibit a considerable uncertainty that is larger than those emerging from re-analysing past events as in the previous section. Because the tsunami hazard is strongly linked to the landslide dynamics, it is important that these uncertainties are reflected in a probabilistic study to weight them in a systematic way. Below, we demonstrate how these landslide uncertainties can be incorporated in a tsunami hazard study, using the example case of Lyngen fjord situated in Northern Norway to demonstrate the framework.

Setting

The study is conducted in the Lyngen fjord situated in Northern Norway, which location is shown in Fig. 2. Four unstable rock slope volumes, all located on the eastern side of the fjord, contain a potential rock slide-induced tsunami threat to the coastal communities along the Lyngen fjord. These rock slides are labelled Jettan, Indre Nordnes, Revdalsfjell 2 and Revdalsfjell 1B. Volumes and motion of these slopes are retrieved by combining satellite- and ground-based InSAR (Eriksen et al. 2017) as well as other physical measurements such as borehole monitoring.

Fig. 2
figure 2

The Lyngen fjord, with the location of the four different unstable rock slopes considered in this study, namely Jettan, Indre Nordnes, Revdalsfjell 2 and Revdalsfjell 1B

Full size image

In this study, we limit the analysis to tsunami hazards failure along these four slopes. Hence, we do not consider the possibility of tsunamis due to other landslides occurring on other slopes. Moreover, the estimation of the landslide recurrence (in terms of annual probability of occurrence) is beyond the scope of this paper. However, we note that the annual probability of occurrence is necessary for quantifying the tsunami hazard. Fortunately, estimates of the probability of occurrence are available from a previous study conducted by NVE (Majala et al. 2016). All the different landslide volumes with associated annual probabilities are retrieved from (Majala et al. 2016) are listed in Table 2. These previously estimated return periods are based on expert judgement by the authors of the report (Majala et al. 2016), using a combination of slope movements measured on the actual unstable slope, combined with information of historical landslide events in the region proximal to the potential landslides (Hermanns et al. 2013). Consequently, areas with larger historical activity and larger slope movements are considered more probable than slopes with lower movement and less evidence for historical events. Table 2 also contains rough estimates of the uncertainty of the annual probability of occurrence, which should cover 90 % of the uncertainty range.

Table 2 Unstable rock slopes with estimated volumes, annual frequencies and uncertainty range in annual frequency for the Lyngen fjord study. The upper and lower ranges in annual frequencies are expected to cover 90% of the uncertainty range. All values are derived from (Majala et al. 2016), except the mean value, which is derived assuming that the annual frequency follows a lognormal distribution
Full size table

Probabilistic and computational framework

Computational framework

A computational framework capable of simulating thousands of landslide tsunami scenario propagation and inundation models is implemented as the basis for the probabilistic analysis. The finite difference tsunami model GloBouss (Løvholt et al. 2008) is used to simulate the open tsunami propagation in the fjord. GloBouss is a Boussinesq-type model and incorporates first-order non-linear terms and optimised frequency dispersion. The tsunami simulations in the fjord are carried out on a local bathymetry with an original grid resolution of 25 m. To find a sufficiently accurate resolution for the tsunami propagation grids, we have performed grid refinement convergence tests. At locations with typical depths similar to the depth at the fjord side of the inundation grids (where the waves from the propagation model are fed into the inundation model), the height of the leading wave for 100-m resolution differs with less than 0.1% compared with 50-m grid resolution. For the second wave, we find a difference of 5–6%. In the present study, the leading wave is generally always the highest one and is the main contributor to the maximum inundation heights. The convergence tests therefore show that for tsunami propagation a grid resolution of 100 m is sufficient, which is hence applied in the simulations. The present simulations are carried out in linear dispersive mode.

Tsunamis are generated by landslide block sources, the methodology is described in, e.g. (Løvholt et al. 2015; Harbitz et al. 2014b; Løvholt et al. 2005) and we only review a few important details here. The landslide has a shape of a rounded rectangular block, having a frontal area Aj, a width Wj and a length Lj, p, where j and p indicate event realisation numbers explained below. The landslide moves along a straight line with a prescribed motion following a sine function in time. The landslide impacts the water with a velocity Um, retarding its motion until it has reached a velocity Rn, where m and n are event realisation numbers. The tsunami generation is then taken into account as a time-dependent volume flux source as the landslide advances as described in (Løvholt et al. 2015).

For estimating possible ranges in landslide velocities and impact velocities as well as frontal areas, we have applied the Voellmy version (Kim 2014) of BingClaw (Kim et al. 2019). We refer to this model as VoellmyClaw below. This information is then used as input for setting up the block model kinematics described in the paragraph above. The reason for using a block source in the probabilistic mass production of simulations instead of using the VoellmyClaw model directly is that we have better control of the shape and kinematic variability of the landslide and because the VoellmyClaw model demands more computational resources and data handling than a block source model. In VoellmyClaw, the main parameters are friction angle and friction from turbulence and viscous drag.

Inundation simulations are also carried out on a series of different coastal locations. The inundation simulations are performed on refined topographic grids having a resolution of 10 m on the finest grid (nested grids). The convergence tests of grids with 5- and 10-m resolution show less than 0.1 % difference of the leading wave in the sea for depth larger than, e.g. 50 m. In the inundated areas, the differences are less than 3 %. From this, we conclude that the 10-m grid gives sufficient accuracy. The MOST model (Titov and Synolakis 1997) is used for these inundation simulations, using a nesting procedure to couple the tsunami propagation model dynamically to the inundation simulation (Løvholt et al. 2010).

Event-tree analysis

An event-tree analysis is set up to analyse tsunami probabilities. The event tree used in the present study is displayed in Fig. 3. The event-tree structure contains branches representing individual conditional landslide event probabilities, which are independent. The full set of branches contains all the landslide parameters considered in the uncertainty analysis. Furthermore, a single branch set contains a complete set of conditional probabilities for a given landslide parameter, for instance, a given landslide volume (indicated as Volr in Fig. 3) or a landslide impact velocity (indicated as Um in Fig. 3).

Fig. 3
figure 3

General structure of the event tree including all branches for the Lyngen study

Full size image

In this study, the present event-tree structure is implemented for estimating the uncertainty in tsunami probability at a given coastal location. The explanation below follows the event tree given Fig. 3 from left to right, hence starting from the top of the tree structure. Except for the probability of the landslide volume, which is given in terms of the annual probability, all conditional probabilities are in the [0 − 1] range, and the sum of all the conditional probabilities for a given parameter is 1. For this, we ended up with the following probability branches:

  • Exceedance probability for inundation height probabilities for a given landslide number q. To retrieve the overall hazard at a given coastal site, the probabilities due to tsunami from all landslides Nq = 4 are combined.

  • The rate λ = P(Volr| landslideq) of a landslide volume failing on slope q. This probability is given in terms of the temporal frequency information.

  • The conditional probability of alternative tsunami generation mechanisms PM = P(Mp| Volr). In this study, this parameter concerns two different alternative models for the tsunami generation. The first alternative is a finite block that slides into the basin from the land. The second alternative is represented by a block with “infinite” length, meaning that the submerged slide length will always be larger than submarine landslide run-out distance. These alternative models for the tsunami generation were introduced due to the imperfection of the numerical model, which is linear in the tsunami generation area. The finite block length captures the landslide volume correctly, but because the rear part generates a large wave through when the rear part of the landslide enters the water, it is not considered fully realistic. Vice versa, the infinite block represents a more realistic generation mechanism, but the volume will be too large.

  • The conditional probability PA = P(Aj| Mp) for a landslide frontal area Aj, given the tsunami generation mechanism Mp (finite or infinite block length).

  • The conditional probability PU = P(Um| Aj) for a landslide impact speed Um, given the frontal area Aj.

  • The conditional probability PR = P(Rn| Um) for the submerged landslide run-out distance Rn, given the landslide speed Um.

The implementation of the event tree assumes that all parameter probabilities are combined, and that all parameter combinations are independent of each other. For instance, this means that all realisations of a landslide velocity Um are combined with all frontal areas Aj, all frontal areas Aj are combined with all landslide volumes Volr and so on. This is clearly a simplification, as the landslide parameters in most cases are dependent on each other. One exception from this assumption is however implemented, namely a dependency between the landslide velocity Um and the run-out distance Rn, where we allow only the highest velocity landslides to reach the largest run-out distances and vice versa. In this case, the full set of branches are not used, and the conditional probabilities are re-normalised by increasing the probabilities of the remaining branch in such a way that the sum of all active realisations is still unity.

We assume that the probability of a single event follows a Poissonian process. First, we define the probability of a single landslide source in the event tree given as the product of all conditional probabilities over the event-tree branch multiplied with the rate λ:

$$ {P}_{q,r,p,j,m,n}=1-{e}^{-\lambda \cdotp {P}_M\cdotp {P}_U\cdotp {P}_A\cdotp {P}_R} $$
(1)

While the overall probability summed over all simulations is given by

$$ P{\left(x,y\right)}_I=1-\prod \limits_{q,r,p,j,m,n}\left(1-\left[{P}_{q,r,p,j,m,n}\cdotp P{\left(x,y\right)}_{q,r,p,j,m,n}\right]\right) $$
(2)

where the exceedance probability of a coastal point being inundated is found from the inundation simulations by taking into account that P(x, y)q, r, p, j, m, n is the binary probability of a point being inundated, which is simply 1 if a point in a scenario simulation is wet, and 0 if the point is dry.

Probabilistic example studies for Lyngen fjord

Discretizing landslide parameters and filling the probabilities

The landslide parameters used in the PTHA event-tree analysis and the related probabilities are here defined by using available knowledge, statistics and experience gained from previous simulations of subaerial landslide tsunami studies. The probabilities were assigned using expert judgement and hence set subjectively through joint discussions among the authors of this paper, based on available data and simulations. Yet, as elaborated below, the chosen parameter values and associated probabilities lend themselves heavily on data where available.

Due to the number of branches in the event tree and because each event requires a significant amount of computational resources, just a few realisations of each landslide parameter is possible for this study. For the submarine run-out distance Rn and the frontal area Aj, we assign five different alternative parameters with five associated probabilities, while for the impact velocity Um we have used three. The parameter values and corresponding probability for these three parameters are found in Table 3.

Table 3 Parameters and corresponding probabilities (P) used in the calculations for all four rock slide locations for the Lyngen study. R is the submarine runout (five values), A the frontal area (five values) and U the impact velocity (three values). Combined with the possibility of two alternative generation mechanisms, this gives 150 different combinations (5 × 5 × 3 × 2) for each location and each volume.
Full size table

The probabilities for the submarine run-out is based on statistics on the available data on the ratio between the fall height (H) and run-out length (L) as a function of the volume for large rockslides in Norway, see (Romstad et al. 2009) and references therein. The five parameter values for the submarine run-out distances are found by taking the mean value, as well as the upper and lower 95 and 67 percentiles from Fig. 4.

Fig. 4
figure 4

The regression analysis of the run-out statistics of rock slides in Norway as explained in the text. H is the fall height, L is the total horizontal run out. The vertical lines A, B and C indicates the volumes 0.8, 5 and 6 Mm3

Full size image

The three different landslide velocities used in the probabilistic analysis are typically lower and upper limits of impact velocities (U1 and U3) for rockslides in Norwegian fjords (Harbitz et al. 2014b), while the mean impact velocity (U2) is calculated from a simple numerical scheme for block slide dynamics labelled the energy line approach.

The process of determining the landslide frontal area and run-out distance are carried out using the following steps: First, the VoellmyClaw model parameters are tuned to reach the mean run-out value (R3 in Table 3) and the velocities from the energy line approach (see above). A special case is Revdalsfjell 1B, which in the case of mean runout do not reach the fjord. In this case, we apply R1 as reference. Area A2 is defined as the total affected area over the coastal cross-section integrated overall simulation times, while A3 is defined as the largest instantaneous area crossing the coastline at any given time. The other values of the front area A1, A4 and A5 are defined as A1=1.5 ∗ A2, A4 = 0.6 ∗ A2 and A5 = 0.4 ∗ A2. The ratio between the height and width of the slide is determined by extracting similar ratios from the VoellmyClaw simulations.

For the probability PM related to cutting the landslide end during impact, we use PM1 = 0.3 (a finite length block) and PM2 = 0.7 (an infinite length block). The reason for the emphasis on the infinitely long block is that the calibration of the landslide dynamics (velocity and runout) in the past tsunami analysis has assumed this approach for fitting previous tsunami simulations towards observations, e.g. (Harbitz et al. 2014b).

We stress again that the annual frequency of each of the different rock slide events was provided from an external study. Hence, establishing the source annual frequencies were not a part of this study. However, the frequencies, as presented in Table 2, were given as median values with associated uncertainty ranges. To derive the mean landslide frequency, which was used as rates in the event tree, we assumed lognormal distribution. Derived mean values are given in Table 2.

Inundation height uncertainty due to single volume landslides

Here, we compute the hazard using all the four source regions listed in Table 2, aggregated over all the landslide parameters and probabilities listed in Table 3. Because we only investigate single volume landslides, we refer to this work as semi-probabilistic (the full probabilistic analysis is carried out in the next section). Figure 5 shows a tsunami hazard for one of the inundation sites (Lyngseidet), for 1/1000 year−1 and 1/5000 year−1 exceedance probabilities. In addition, we show the maximum inundated area aggregated for all simulations. Due to the steep slopes in the Lyngen fjord, the contour maps for the 1/1000 year−1 and 1/5000 year−1 exceedance probabilities coincide in many places. This is due to the steep slope controlling the inundation distance, which does not exclude differences in inundation heights. On the other hand, we see that the most extreme scenarios in the probabilistic analysis by far exceed the inundation limits for these return periods.

Fig. 5
figure 5

Inundation height exceedance probabilities for the Lyngseidet site using the fixed volume approach. Three different levels are considered: 1/1000 year−1 exceedance probability (red line), 1/5000 year−1 exceedance probability (black line) and the maximum inundated area overall simulations (green line)

Full size image

In this investigation, a total of 31 probabilistic high-resolution inundation maps were computed, often with the same characteristic behaviour as revealed for Lyngseidet. However, only a few sites like Lyngseidet were located in the immediate vicinity of the identified landslide prone sites shown in Fig. 2. Hence, many locations faced relatively smaller hazard. To this end, we aggregated the hazard for all the sites into joint hazard maps. These hazard maps are shown in Figs. 6 and 7 for 1/1000 year−1 and 1/5000 year−1 exceedance probabilities, respectively. As shown for the 1/1000 year−1 exceedance probability, tsunami heights in the vicinity of the landslide areas (along a stretch of 5–10 km north and south from two northernmost landslides, say), we find maximum inundation heights in the order of 4–8 m. Further afield, the inundation heights rapidly drop off to 3 m and less. The 1/5000 year−1 exceedance probability inundation heights are considerably larger, with associated inundation heights up to 10 m.

Fig. 6
figure 6

1/1000 year−1 exceedance probability for maximum inundation heights aggregated for all inundation sites using the fixed volume approach. The yellow markers indicate the location of the four rock slides

Full size image
Fig. 7
figure 7

1/5000 year−1 exceedance probability for maximum inundation heights aggregated for all inundation sites using the fixed volume approach

Full size image

The 1/5000 year−1 hazard map corresponds to the strictest of the two hazard maps applied to Norwegian building legislation. Previously, the authors of this paper have carried out also deterministic tsunami hazard analysis (unpublished). When comparing the new semi-probabilistic 1/5000 year−1 hazard results with the deterministic ones, we typically find that the deterministic results provide at least 20–30% larger inundation heights than the semi-probabilistic method (results not shown). Hence, by introducing the probabilistic method to weight different alternative models, there is a tendency in this case that less conservative assumptions are made. Consequently, areas previously located in hazard zones could be freed.

We notice that the hazard is considerably larger in sites close to the Jettan landslide, which is considered having the highest probability to fail compared with the three other landslides, as revealed by Table 2. Naturally, the hazard levels are sensitive to the mean landslide recurrence rates because the landslide volumes are fairly similar. We emphasise again that these mean landslide rates have been determined expert judgement by a group of engineering geologists carefully assessing each slope (Majala et al. 2016). In the present semi-probabilistic approach, the uncertainty related to the return period and the volume have not been quantified, as only the mean value of the return period was used. In the more general PTHA carried out in the next section, we explore briefly the sensitivity to these assumptions.

Synthetic PTHA study containing multiple volumes

Contrary to the analysis in the previous section, we will in the following adopt a PTHA study where multiple landslide volumes contribute to the hazard. We consider only inundation to the Lyngseidet site, already analysed using the semi-probabilistic approach (see Fig. 5). As the inundation at this site is heavily dominated by tsunamis caused by Jettan and Indre Nordnes slope failures, we further consider only these two source regions. Moreover, we assume exactly the same landslide parameters and related parameters as we did for the volumes considered in the previous section. We realise that this assumption would likely be violated in nature, as the landslide velocities, frontal areas and run-out distances from smaller landslides are likely to be smaller than for the large landslides. In any case, we adopt this assumption because it simplifies our analysis, but also because the main purpose here is to demonstrate the method.

In the following, we use the tapered Pareto distribution (see, e.g. Geist and Parsons, 2014) to generate magnitude frequency distributions (MFDs) relating the volumes to the probability of occurrence for the Jettan and Indre Nordnes sites. The tapered Pareto distribution used here is given by

$$ P\left(V\le v\right)=P0{\left(\frac{V}{V_{\mathrm{max}}}\right)}^{\beta}\cdot {e}^{\frac{V-{V}_{\mathrm{max}}}{V_c}} $$
(3)

where P0 is a normalisation factor, V the landslide volume, Vmax the largest volume, Vc the value for the onset of the tapering and β as steepness factor. P0 is adjusted to make the integrated annual probability for all volumes in the analysis being equal to 1/50. The 1/50 factor was chosen subjectively but was necessary for demonstrating the methodology. We remark that this value was kept low to comply with the lack of historical records of large rock slides. A value of 1/50 seems plausible due to the lack of large historical events in the areas and could even have been given a lower value on this basis.

The parameters used in the present MFDs are listed in Table 4. These parameters make use of three different MFDs for each of the two landslide areas. In all cases, we adjust the parameter β to obtain the desired rates for the most voluminous landslides. MFD1 gives similar rates for the largest volumes (V ≥ 6⋅106 Mm3) as the single-volume analysis in the previous section. MFD2 gives a larger probability to the landslides with V≥ 6⋅106 Mm3 and a lower probability to the smaller volumes. The opposite is the case for MFD3. To establish landslide scenarios, we discretise each MFD into eight different volumes separated logarithmically. Then, we simulate the landslide dynamics and tsunami generation as described in the “Probabilistic and computational framework” section using the landslide volumes identified for the MFDs, and all the different scenario parameterisations given in Table 3. All MFDs used are depicted in Fig. 8, with the scenario volumes given on the x-axis and the rate on the y-axis.

Table 4 Values used for creating synthetic MFDs. The total annual rates for volumes exceeding 6 Mm3 are also listed
Full size table
Fig. 8
figure 8

Synthetic magnitude frequency distribution curves for landslide volumes for the Jettan landslide (upper panel) and the Indre Nordnes landslide (lower panel). The different MFDs are used in PTHA to illuminate the sensitivity to modelled annual source rates. The markers in the figures show the rates used for the individual scenario realisations

Full size image

Results using the PTHA approach are summarised in Fig. 9 for two exceedance probabilities, namely 1/1000 and 1/5000 year−1. The results are also compared with the tsunami inundation maps using the semi-probabilistic analysis represented by a single volume landslide in each location from the “Inundation height uncertainty due to single volume landslides” section. From the results, we make the following observations, recalling that MFD1 puts similar weight on the largest volumes as in the semi-probabilistic analysis, that MFD2 puts more weight on the largest volumes and MFD3 puts more weight on the smaller volumes: For the 1/1000 year−1 exceedance probability map, we see clearly that results using MFD2 cover a larger area than the other analyses. Furthermore, it is observed that MFD1 cover a slightly larger area than the semi-probabilistic analysis, while the MFD3-based analysis gives almost the same inundated area as the semi-probabilistic analysis, albeit with slightly smaller inundation distance. When we move to comparing the 1/5000 year−1 exceedance probability map, we find less pronounced differences between the four different analyses, with the exception that the MFD2 analysis provides a larger inundated area than the three other analysis sets (which are almost identical).

Fig. 9
figure 9

Comparing tsunami inundation maps using different synthetic MFDs. Left panel 1/1000 year−1 exceedance probability. Right panel 1/5000 year−1 exceedance probability. In both cases, we compare the hazard maps with results obtained using the fixed volume approach

Full size image

It is clear that the results are most sensitive to the choice of MFD for the 1/1000 year−1 exceedance probability compared with the 1/5000 year−1 one. Comparing the results using different MFDs with the semi-probabilistic analysis did not produce dramatically different results in the present case. However, part of the explanation for the small differences also lies in the strong topography control on the probability contours in steep topography. In a setting with a flatter topography, the differences could possibly be more pronounced. Moreover, the uncertainty range we explore in the different synthetic MFDs is not very large given the large lack of knowledge on the failure probabilities. As the uncertainty is directly linked to MFD, exploring larger uncertainties would naturally also be reflected in the hazard maps. While it is beyond the scope of this paper to analyse the probability of source recurrence, it is obvious that this parameter has a first-order effect on the hazard.

In the present framework, we carried out probabilistic tsunami hazard analysis using single MFD curves as input. A more general approach would be to treat such models as alternative realisations, represent the hazard with related uncertainty limits. Such an approach is already implemented for earthquake tsunamis (Selva et al. 2016), and it would in principle be possible also to apply it to landslide PTHA in the future.

Concluding remarks

In this paper, we have reviewed how landslide and tsunami data can be used to better constrain landslide dynamics uncertainty, and in turn, landslide tsunami probability. Past landslides can be used to train models to provide bounded envelopes for landslide parameters. In cases where only landslide run-out data are available, these uncertain parameter bounds are usually rather broad and sometimes not even properly constrained. It is further shown that tsunami data can dramatically reduce the interpreted uncertainty bounds for the landslide parameters. The reason for this is that the tsunami is an indirect footprint of the landslide motion.

For prognostic analysis of landslide tsunamis, we need to resort to both data as well as to subjective choices for quantifying the landslide tsunami hazard. In this paper, a framework for aggregating probabilities from different uncertain parameters is combined using a simple event-tree analysis. In this analysis, probabilities stemming from both empirical data and judgement are combined in the same framework. A more rigorous framework, separating treatment of aleatory (inherent) uncertainty and epistemic (lack of knowledge) uncertainty is not attempted here. However, we admit that more rigorous uncertainty treatment should be attempted for PTHA from landslide sources in the future.

An example of a LPTHA application for the Lyngen fjord area in Northern Norway is presented and discussed. The results from this LPTHA is largely controlled by the annual rates of occurrence from the subaerial landslides that are input to the analysis. However, there is also large variability between the different landslide parameter realisations, which gives rise to significant differences in the hazard. To this end, the worst-case scenario within the set of events comprising the PTHA study provided a much larger inundation area compared with the 1/1000 year−1 and 5000 year−1 year exceedance probability subject to Norwegian legislation. In general, it was found that by using a probabilistic analysis rather than a deterministic one, our LPTHA analysis tended to be less conservative. Hence, new inundation limits were less strict than previously enforced by the legislation. Another finding from the probabilistic analysis was that the largest landslide volume, at the corner frequency of the magnitude frequency distribution, say, largely determined the tsunami hazard at the return periods of interest. This deserves to be better investigated though.

In summary, the probabilistic landslide tsunami hazard analysis is still in its infancy. However, advancements are being made in better understanding landslide tsunami uncertainty. The probabilistic framework is in place to embed such uncertainties into LPTHA, but we still need to use judgement for estimating probabilities. By better making use of past landslide and tsunami data, as well as high-performance computing resources, we expect LPTHA to provide more mature input to future hazard assessments. On the other hand, landslide tsunamis are infrequent events, and the hazard should always be expected to carry large uncertainties.