Introduction

On January 29, 2014, a submarine landslide and a resulting tsunami caused damage to port facilities and nearshore infrastructure at Statland, county of Nord-Trøndelag, Norway. The Statland event is the last one out of a series of tsunamis that have inundated Norwegian coastlines. These tsunamis are generated by landslides, although earthquakes may, in certain cases, also act as the trigger. Through history, the most damaging tsunamis in Norway are those induced by subaerial rockslides impacting fjords or lakes. Three subaerial rockslides in Tafjord 1934 and Loen 1905 and 1936 caused altogether 174 fatalities in the last century (Harbitz et al. 1993). The second example of tsunamigenic landslides is submarine clay-rich landslides that commonly involve a retrogressive landslide development (Løvholt et al. 2016). Perhaps the most well-known example is the 8150-BP Storegga Slide tsunami that involved transoceanic propagation inducing run-up heights that exceeded 5 m in Scotland and possibly in Denmark (e.g., Smith et al. 2004; Fruergaard et al. 2015) and 10 m in Norway, Shetland, and the Faroe Islands (e.g., Bondevik et al. 2005).

While the Storegga Slide involved an enormous volume, several smaller cases of submarine, multistage landslides and associated tsunamis in nearshore environments such as fjords, lakes, and estuaries have previously been reported in Norway. Examples of such tsunamis are the ones in Trondheim harbor (1888), Orkdalsfjorden (1930), Sokkelvik (1959), Rissa (1978), and Balsfjord (1988); see L’Heureux et al. (2013, 2014) and references therein. It has further been speculated that a submarine landslide caused tsunami inundation that damaged Stone Age settlements in Rennesøy, SW Norway (Bøe et al. 2007). The Statland event falls into this category of nearshore submarine landslide tsunamis.

It is vital to increase our understanding of the processes involved in landslide failure, its dynamics, and tsunami generation in order to ensure a safer and more reliable urban development in nearshore environments. However, few well-documented cases of landslide tsunamis exist. Some of the best-studied examples are the fully submerged 1998 Papua New Guinea tsunami (Kawata et al. 1999; Synolakis et al. 2002; Tappin et al. 2008) and the 1929 Grand Banks tsunami (Fine et al. 2005). Other well-known examples involving smaller volumes and coastal failure are the 1979 Nice event (Assier-Rzadkiewicz et al. 2000), the twin tsunami in Haiti (Fritz et al. 2013), the 1994 Skagway tsunami (Kulikov et al. 1996), and a series of landslide tsunamis following the 1964 Great Alaska earthquake (Parsons et al. 2014). The latter events are similar to the Statland one in that they involved failure of soft marine sediments close to the shoreline. As will be demonstrated, the availability of both pre- and post-landslide seafloor data as well as inundation measurements make the Statland event a good test case toward a better understanding of landslide-tsunami processes. The work presented in this paper was part of an investigation to find the causes of the Statland landslide (NVE 2014).

The deltaic deposits outside the village of Statland have accumulated over the years along the margin of Namsfjorden. The deposits consist mainly of loose sands and silts lying on marine and partly sensitive clays (quick clays). Swath bathymetry data were collected by Seascan by means of a 250-kHz GeoSwath system prior to and after the landslide. The data show that up to 20 m of sediments over an area of 21,000 m2 were displaced due to the landslide (NVE 2014). The total volume of sediments amounts to 350,000–400,000 m3. Seafloor mapping and seismic profiling performed by the Geological Survey of Norway (NGU) further show that the run-out distance of the landslide debris was 1300 m in Namsfjorden, the landslide drop height about 380 m, and the deposit thickness 1–3 m (NVE 2014).

In the present paper, we describe eyewitness observations as well as results of the field survey that was conducted 2 days after the landslide. Secondly, we describe the landslide and tsunami modeling for two different scenarios. Finally, the numerical results from the tsunami simulations are compared with eyewitness accounts and field survey measurements to make inferences on the evolution of the landslide dynamics.

Field survey and eyewitness observations

The Statland landslide occurred at 04:31 p.m., January 29, 2014. The mean high tide in this area is about 0.8 m above mean sea level (MSL). Fortunately, the water level was almost at the lowest astronomic tide (LAT), about 1.4 m below MSL at the time of the tsunami impact. A field survey covering the tsunami inundation as well as a qualitative description of the damage due to the tsunami were conducted by the Norwegian Geotechnical Institute (NGI) on January 31, 2 days after the landslide occurred. The field investigation revealed that the measured maximum run-up height was 8.2 m above MSL (9.6 m above contemporary sea level) in the Sagvika bay; see location 2 in Fig. 1.

Fig. 1
figure 1

Overview over Statland. The initial landslide is marked with red color (labeled 1). 2 The highest run-up, about 10 m above the contemporary sea level. 3 Several boathouses destroyed at the headland Kvalvikskjæret. 4 and 5 Water reaching a few meters inland (flow depth at point 5 probably only a few tens of cm). The letters A–H refer to the photos in Fig. 2. The location of Statland is shown with a red bullet in the inlet map of Norway

Full size image

The tsunami involved no fatalities, although one person got minor injuries. The injured person was inside the large industrial building located in the middle of Sagvika when it was hit by water from the tsunami flowing 1 m deep. As a consequence, he fell and was pushed 7–8 m along the floor. The industrial building was severely damaged, while 12 boathouses were totally destroyed (see photos from the field survey in Fig. 2). In addition, the water supply to the land-based aquaculture broke (located east of point 5 in Fig. 1), and a section of the road close to the landslide area was swept into the fjord. Further, dockages outside four boathouses were destroyed after being hit by the tsunami (west of point 5), and the wooden dock located centrally in Sagvika was crushed and the floating docks further north in Sagvika were washed on land. Several boats (cabin cruisers or smaller) sank both in Sagvika and Djupvika.

Fig. 2
figure 2

Photos from the field survey. a The landslide scar. b The industrial building in Sagvika with the floating dock washed on land on the opposite side of the bay Sagvika in front. c The road behind the industry building in Sagvika. d Boathouses close to the industrial building in Sagvika. e Outside the landslide scar, between Sagvika and Djupvika. f Boathouses at Kvalvikskjæret. g Destroyed dockings in Djupvika. h Run-up in Djupvika. See Fig. 1 for locations. In the photos, the water level is about 2.1–2.4 m above contemporary sea level at the time the landslide was triggered (or about 0.7 to 1.0 m above MSL)

Full size image

In the 1st floor of the local grocery (located at the headland northeast of point 2 in Fig. 1), a woman witnessed water flowing along all four walls outside of the building, but due to the limited flow depth, significant structural damage was not sustained. Several eyewitnesses observed waves entering Djupvika bay. First, a 3–4-m-high wave propagated toward the bay. The first positive wave was followed by a larger depression, reported by eyewitnesses to be 4–6 m below MSL 1 min after the leading wave peak entered the bay. A minute later, the water started to refill the bay with water. A dramatic and huge sound due to the water motion was also reported, and the entire bay was filled with large vortices.

The tsunami trimline in Sagvika and further north was traced successfully with particularly dense sampling along the roads in the area (e.g., photo C in Fig. 2). However, in the eastern part of the survey area, the run-up height was lower, and due to the low contemporary water level the run-up was limited to the swash zone. Hence, onshore traces were not found (e.g., photo E in Fig. 2) except for two smaller areas at point 4 and point 5 (Fig. 1). At Kvalvikskjæret (point 3, Fig. 1), some of the destroyed boathouses were swept over the headland and into the sea on the northern side.

The measured run-up heights are shown in Fig. 3. In Sagvika, the maximum run-up (measured from MSL) was 8.2 m, at Kvalvikskjæret 4–7 m, and in Djupvika 2.5–3.5 m.

Fig. 3
figure 3

Measured run-up heights along the trimline using a GPS and the heights from laser scanning along the track. Height is given in meters and is measured above the vertical datum NN1954. NN1954 at Statland is 17 cm above MSL. The contemporary water level was 1.4 m below NN1954. The landslide is marked with the dashed black line. The coordinates are given in UTM zone 32

Full size image

Landslide dynamics and run-out

We used the BING model (Imran et al. 2001) for modeling the landslide dynamics. The BING model assumes depth-averaged visco-plastic Bingham rheology, where the landslide evolves as a combined plug flow riding on top of a shear layer. The plug and shear layer thicknesses are determined by the slope and the yield stress of the material. The model is quasi two-dimensional (2D) and assumes that the sediments disintegrate and liquefy instantly upon failure. The model takes the buoyancy into account, but not the added mass and viscous drag. Omitting the drag and added mass means higher velocities and accelerations (De Blasio et al. 2003) in the BING model, which may in turn influence the wave generation unless it is accounted for. More details about the BING model, including the governing equations and derivations, are included in Imran et al. (2001). The BING model has also been used in other studies to describe submarine landslide dynamics in tsunami modeling (e.g., L’Heureux et al. 2012; Glimsdal et al. 2013; Harbitz et al. 2014; Løvholt et al. 2014).

We propose two scenarios representing different hypotheses for how the Statland landslide failed and evolved. The first and simplest scenario assumes that the entire volume is released in one event (see Fig. 4). However, morphological analysis shows that the landslide most likely developed in several stages (see Fig. 5). The second scenario therefore assumes a two-stage evolution of the Statland landslide. The landslide dimensions for the different scenarios are given in Table 1 while material parameters for BING and main model outputs are found in Table 2. The input parameters to BING are based on empirical correlations with available geotechnical parameters (see Locat and Demers 1988). The resulting yield strength in BING was found to be similar to that of an average remolded undrained shear strength for the failed sediment volume in each scenario.

Fig. 4
figure 4

The initial thickness of the Statland landslide before release based on the difference between the swath bathymetry data collected prior to and after the landslide. The sediment thickness is given in meters. The coordinates are given in UTM zone 32

Full size image
Fig. 5
figure 5

The suggested two-stage evolution of the Statland landslide. The initial landslide is released first (stage 1) and is followed by the larger quick clay landslide (stage 2). The bathymetry shown in the figure is from the seabed scanning after the landslide. The blue line represent the shoreline prior to the landslide

Full size image
Table 1 Landslide parameters for the different stages. Only the maximum landslide thickness and the initial landslide length are used as input parameters to the BING model
Full size table
Table 2 Input parameters and results from the BING model for the three scenarios
Full size table

For the first scenario, where the entire volume is released simultaneously, the modeled thickness of the landslide deposits was found to be too high compared to the observations (1–10 m modeled versus 1–3 m found from the seafloor mapping). In the second scenario, a smaller frontal part of the sediments was released (labeled “initial slide”) before the rest of the landslide consisting of quick clay was set in motion (labeled “quick-clay slide”). The seafloor morphology reveals that the initial landslide debris first traveled north before it continued its descent to Namsfjorden to the northeast. The modeled landslide deposit for the second scenario is in the range of 1–2 m, which agrees well with results from seismic reflection surveys (NVE 2014). We primarily focus on the two-stage scenario below.

A fixed shape box moving with pre-scribed kinematics was used as tsunami sources in the numerical tsunami model. Løvholt et al. (2015) exemplified that late post-failure deformation has limited effect on the tsunami generation. With the given uncertainties in question regarding staged failure and material parameters, a fixed shaped landslide was therefore considered sufficient in the present example. We further assume a sinusoidal prescribed analytical velocity profile for the fixed shaped box (see Fig. 6; the procedure is adopted from Løvholt et al. 2005). The landslide shape and kinematics for the box-shaped tsunami sources are adopted from slide velocities from the BING model. For the prescribed velocity profile, we apply slightly lower initial accelerations compared to the BING runs, to account for the lack of added mass in BING (omitting added mass typically overestimates the acceleration by 20–30 %; see Watts 2000). Similarly, due to the lack of hydrodynamic drag in the BING model, the maximum velocity in the applied velocity profile for the tsunami model is set to be 20 % less than the one obtained directly from BING. The choice of 20 % reduction is based on our general experience for a landslide such as Statland with relatively short run-out. The 20 % reduction in landslide speed was confirmed by idealized sensitivity tests using a simple box-slide model with skin drag (results not shown).

Fig. 6
figure 6

Landslide velocity. Labels “bing” and “mod” are the results from the BING model and modified results from BING used in the tsunami modeling, respectively. Stage 1 (initial landslide) and stage 2 (quick clay landslide) refer to the scenario with a two-stage release of the landslide as explained in the text

Full size image

The second scenario consists of an initial landslide (stage 1), which is bell-shaped with a maximum height of 15 m, and the quick clay slide (stage 2) which is modeled as a flexible rounded box. The trajectories for the two slide phases differ, as stage 1 follows a curved trajectory (Fig. 7), while stage 2 follows a straight line (Fig. 8). Since the largest wave impact at Sagvika and the strongest withdrawal of water in Djupvika are mainly generated in shallow water, we are led to believe that the difference in slide paths for stages 1 and 2 is not of primary importance. In the present scenario, the quick clay landslide is released 10 s after the initial landslide, i.e., when the already advancing initial landslide is turning eastward (see the curved path in Fig. 7). The 10-s time lag is only a qualified guess, but as we discuss below, numerical tests confirm that the time lag has a limited effect on the maximum run-up.

Fig. 7
figure 7

Modeling of the initial landslide (stage 1). The thickness of the initial landslide is given in meters. The landslide path is marked with a black line, while the boundary of the entire landslide before release is drawn with a red line. The coordinates are given in UTM zone 32

Full size image
Fig. 8
figure 8

Input from landslide modeling to the tsunami propagation model (vertical flux per unit area due to landslide in m/s). Maximum positive value (front of initial landslide) is about 16 m/s. The timing of the snapshots are (from upper left to lower right) 6, 15, 25, and 80 s after release of the initial landslide. The path for the initial landslide and the quick clay landslide is marked with a dotted red line and a solid black line, respectively, while the boundary of the entire landslide before release is drawn with a red line

Full size image

The tsunami propagation model uses the rate of change in landslide thickness, representing a change in volume flux at each time step, to couple the landslide dynamics to the tsunami generation (see Løvholt et al. 2015). Examples of snapshots of the volume flux at four different times are given in Fig. 8, indirectly indicating the position of the front (positive flux) and tail (negative flux) of the landslide. The maximum flux per unit area for the uplift is about 16 m/s.

Numerical tests for the two-stage scenario show that the main contribution to the maximum run-up in Sagvika originates from the front of the initial landslide (stage 1). Correspondingly, the main contribution to the large withdrawal of water in Djupvika originates from the tail of the quick clay landslide. We therefore simplify our modeling of the generation, by assuming that the positive wave originates only from front stage 1 and that the drawdown due originates only from stage 2. These two waves are then superimposed to model the generation due to both slide phases. The full landslide volume after release is bell-shaped in the front with gradually decaying tail (results not shown). The tail (the quick clay landslide) is initially about 250 m long. Quick clay landslides are often released in a retrogressive style with time lags between each release. To model the large wave trough which complies with eyewitness observations in Djupvika, it is clear that the retrogressive back stepping was relatively fast. To incorporate a rapid landslide retrogressive mass release in a first approximation for the wave generation, we model the quick clay landslide as a flexible rounded box with gradually decaying tail.

Tsunami modeling

For modeling the deep-water tsunami generation and propagation, we have applied the GloBouss model (Løvholt et al. 2008; Pedersen and Løvholt 2008). This model includes non-linearities as well as higher-order frequency dispersion, but not the dry land inundation. In this study, the GloBouss model is run with the simple linear water approximation for the deep water propagation, as dispersion and non-linearity were found mostly negligible. We use a prescribed rounded box approximation for including the landslide as a source in the tsunami simulations as outlined in the previous section. For nearshore propagation and inundation, we apply the MOST model (Titov and Synolakis 1995, 1998; Titov and Gonzalez 1997), taking non-linearities in the shallow water propagation (including potential wave breaking) into account. By using a one-way nesting procedure, MOST reads the surface elevation and velocity components from the propagation model (GloBouss) over the model boundaries at each time step. This enables a swift one-way nesting of wave propagation with the run-up model (Løvholt et al. 2010).

The potential for tsunami generation by submarine landslides is highly dependent on the Froude number (landslide velocity-to-wave speed ratio; for a discussion, see Løvholt et al. 2015). A submarine landslide is most tsunamigenic for a Froude number close to unity. The linear hydrostatic wave propagation speed increases with the square root of the water depth. For the Statland landslide, the Froude number is close to unity for the first 400 m of the run-out. In deeper water, the Froude number decays quickly, as the generated waves start to outrun the landslide and the buildup of the waves are strongly reduced. In addition, the initial acceleration may play a role when it comes to the interaction (cancelation) of the frontal and rear waves. However, the high Froude number and short propagation distance mean less interaction between the frontal (positive) and rear (negative) waves. Therefore, the initial acceleration is most likely of less importance here.

In Fig. 9, we show four snapshots of the surface elevation induced by the two-stage landslide using the tsunami propagation model after 10, 20, 30, and 40 s. After 40 s, the wave trough enters the bay Djupvika. Note that the waves are not allowed to inundate dry land in the present model, but is reflected at the shoreline (modeled as a vertical no-flux boundary). The accuracy of these results is confirmed by grid refinement convergence tests for both GloBouss and MOST modeling parts of the current case study (NVE 2014).

Fig. 9
figure 9

Snapshots of the surface elevation for the tsunami propagation model GloBouss after 10 (from top left), 20, 30, and 40 s. The deep wave trough is entering Djupvika at 40 s. The path for the initial landslide is marked with a dotted black line, while the boundary of the entire landslide before release is drawn with a red line. The coordinates are given in UTM zone 32. The boundary of the two areas where the inundation is calculated is marked with black boxes in the upper left panel

Full size image

If the landslide moves into the boundary of the computational domain for the run-up model, the run-up model will not get information about the generated waves inside the domain, because the input from the propagation model is conveyed to the run-up model through the boundaries. For this reason, the domains of the run-up model do not coincide with the landslide area. We have therefore split the run-up simulations into two different domains, one to the northwest and one to the southeast of the landslide (see Fig. 9). For each area, the run-up model calculates the tsunami inundation by using nested grids on three levels, with grid resolutions of 8, 4, and 1 m, respectively. Note that due to short distances between the source area and the coastlines, the coupling between GloBouss and MOST may lead to false waves when waves are reflected after inundation and enter the boundaries/coupling zone between the models. However, the extrema are not influenced by this effect.

Results for the first scenario show that if the entire landslide is released at the same time, it can be inferred that the northward velocity and acceleration will be too slow (in the very first stage in Sagvika), because the maximum run-up in Sagvika is underestimated. However, a one-stage progression mimics the withdrawal of water in Djupvika observed by eyewitnesses. Better match for the maximum run-up in Sagvika is achieved with a two-stage scenario (and at the same time close to identical results for the withdrawal in Djupvika as for the one-stage scenario). In Fig. 10, the measured trimline from the field survey is compared to the trimline obtained from the numerical model using the two-stage scenario. Good correspondence is found in almost all areas, except for the east-most area (where the run-up height is about 20 % too high) and the headland to the north of the bay Sagvika (where the run-up height is about 50 % too high). This substantiates further that a two-stage progression is the most likely process, as also indicated by the seafloor mapping and seismic profiling (NVE 2014). Even though the modeling shows good agreement with the measured tsunami run-up height and the eyewitness observations in Djupvika, the deviations north of Sagvika show that all the details of the landslide progression are not fully captured.

Fig. 10
figure 10

Observed trimline (red) compared to the computed one (black). The boundary of the entire landslide before release is drawn with a dashed black line. The modeled tsunami run-up heights agree with the measured run-up in almost all areas except for the two areas indicated by blue and green arrows. At the green arrow, the modeled run-up is about 20 % too high (about 1 m) while at the blue arrow, the differences are about 50 % (about 4 m). The coordinates are given in UTM zone 32

Full size image

The time lag between the initial and quick clay landslide will not influence the maximum run-up significantly since the main contribution to the maximum run-up is related to the northwest propagation of the initial landslide and not by the quick clay landslide propagating more northeast. This is confirmed by numerical tests (not reported).

Conclusions

The Statland landslide most likely developed in two pronounced separate stages. Both morphological seafloor interpretations and results from combined numerical modeling of the landslide dynamics and tsunami generation and run-up support this conclusion. Analysis of the seafloor morphology indicates that the initial landslide debris must first have traveled northwest before it continued its descent to Namsfjorden to the northeast. The larger (stage 2) quick clay landslide was released after a short time delay and is suggested to have followed a different trajectory more directly into the fjord than the initial part of the landslide. The modeled tsunami run-up heights agree well with the measured run-up and eyewitness observations in almost all areas, except for the eastern-most area as well as at the headland to the north of Sagvika bay. Even though the modeling shows good agreement, the larger deviations north of Sagvika show that some elements of the landslide progression and tsunami generation are not fully captured.