Articles | Volume 22, issue 12
Nat. Hazards Earth Syst. Sci., 22, 3839–3858, 2022

Special issue: Tsunamis: from source processes to coastal hazard and...

Nat. Hazards Earth Syst. Sci., 22, 3839–3858, 2022
Research article
30 Nov 2022
Research article | 30 Nov 2022

Simulation of tsunami induced by a submarine landslide in a glaciomarine margin: the case of Storfjorden LS-1 (southwestern Svalbard Islands)

Simulation of tsunami induced by a submarine landslide in a glaciomarine margin: the case of Storfjorden LS-1 (southwestern Svalbard Islands)
María Teresa Pedrosa-González1, José Manuel González-Vida2, Jesús Galindo-Záldivar1,4, Sergio Ortega2, Manuel Jesús Castro3, David Casas5, and Gemma Ercilla5 María Teresa Pedrosa-González et al.
  • 1Departamento de Geodinámica, Universidad de Granada, 18071 Granada, Spain
  • 2Departamento de Matemática Aplicada, Escuela de Ingenierías Industriales, Universidad de Málaga, 29071 Málaga, Spain
  • 3Departamento de Análisis Matemático, Estadística e Investigación Operativa, y Matemática Aplicada, Facultad de Ciencias, Universidad de Málaga, 29071 Málaga, Spain
  • 4Departamento de Geociencias Marinas, Instituto Andaluz de Ciencias de la Tierra, Consejo Superior de Investigaciones Científicas y Universidad de Granada (CSIC-UGR), 18100 Granada, Spain
  • 5Marine Geosciences, Institut de Ciències del Mar (ICM), Consejo Superior de Investigaciones Cientificas (CSIC), 08003 Barcelona, Spain

Correspondence: Gemma Ercilla (


A modelling approach to understand the tsunamigenic potentiality of submarine landslides will provide new perspectives on tsunami hazard threat, mostly in polar margins where global climatic change and its related ocean warming may induce future landslides. Here, we use the L-ML-HySEA (Landslide Multilayer Hyperbolic Systems and Efficient Algorithms) numerical model, including wave dispersion, to provide new insights into factors controlling the tsunami characteristics triggered by the Storfjorden LS-1 landslide (southwestern Svalbard). Tsunami waves, determined mainly by the sliding mechanism and the bathymetry, consist of two initial wave dipoles, with troughs to the northeast (Spitsbergen and towards the continent) and crests to the south (seawards) and southwest (Bear Island), reaching more than 3 m of amplitude above the landslide and finally merging into a single wave dipole. The tsunami wave propagation and its coastal impact are governed by the Storfjorden and Kveithola glacial troughs and by the bordering Spitsbergen Bank, which shape the continental shelf. This local bathymetry controls the direction of propagation with a crescent shape front, in plan view, and is responsible for shoaling effects of amplitude values (4.2 m in trough to 4.3 m in crest), amplification (3.7 m in trough to 4 m in crest) and diffraction of the tsunami waves, as well as influencing their coastal impact times.

1 Introduction

Submarine landslides represent one of the most common potential offshore geohazards in the continental slopes of the northern high-latitude margins (Elverhøi et al., 2002; Lee et al., 2009). There, the slope failures are essentially focused at their trough mouth fans (Dowdeswell et al., 2008; Rebesco et al., 2014; Llopart et al., 2016; Ercilla et al., 2022, and references therein). Some landslides may also cause tsunamis, as has been evidenced for the Storegga landslide (3000 km3, at 8.1 kyr), with wave amplitudes of  20 m at the shore (Bondevik et al., 2003; Haflidason et al., 2005; Kvalstad et al., 2005), or the Hinlopen landslide (1150 km3, at 30 kyr), with wave amplitudes of  40 m (a.s.l.) (Winkelmann et al., 2008; Vanneste et al., 2006). In the Fram Strait (500 to 1000 km3), located between the main ice retreat areas of Greenland and Svalbard, waves of up to 5.6 m were triggered (Berndt et al., 2009). The factors controlling slope failures in the northern high-latitude margins are still not fully understood. The most common causal factors are the interlayering of underconsolidated glacially derived sediments and low-permeability interglacial hemipelagic clay rich layers, combined with tectonic- and isostatic-related seismicity and/or gas hydrate dissociation (Canals et al., 2004; Kvalstad et al., 2005;Sierro et al., 2009; García et al., 2011; Casas et al., 2013; Vanneste et al., 2014; Moernaut et al., 2017; Llopart et al., 2019).

The understanding of the tsunamigenic potential of submarine landslides still needs to be improved (Chiocci and Ridente, 2011; Løvholt et al., 2020). In this sense back analysis of specific events is commonly used to advance their understanding, as well as to contribute to the hazard assessments of future landslides (Macías et al., 2015; Rodríguez-Morata et al., 2019; Sun and Leslie, 2020). In the European northern high-latitude margins, e.g. Svalbard and Greenland coasts, the tsunami threat has been assessed for a few past landslides. It is important to point out that the tsunami geohazard associated with Holocene landslides such as the Bjørnøyrenna (Laberg and Vorren, 1993), which is the largest landslide (volume of about 1100 km3) in the Barents Sea continental margin; the Nyk landslide in central Norway (Lindberg et al., 2004); or the giant Andøya landslide (northeastern Norwegian–Greenland Sea) (Bugge et al., 1987; Laberg et al., 2000), have not been accurately modelled. The tsunami modelling of recent landslides in this region is important because it will allow us to infer the tsunami potential of future landslides due to climate change and its related ocean warming. Both those interconnected issues are significantly affecting the northern high-latitude margins and may contribute to an increase in the occurrence of submarine landslides, both large and small, in the nearby future, mostly due to gas hydrate dissociation and isostatic-rebound-related earthquakes (Maslin et al., 1998; Tappin, 2010; Urlaub et al., 2013).

Today, the archipelago of Svalbard is one the fastest warming areas of the Arctic Ocean, experiencing an increase in the melting of their glaciers and a rise in the temperature of ocean water circulating along its continental margin (Meleshko et al., 2004; Førland et al., 2013; Skogseth, 2020). This fact may provide adequate conditions to trigger unloading earthquakes and to increase pore water pressure by gas hydrate breakdown, which can destabilize slope sediments (Solheim et al., 2005; Berndt et al., 2009), i.e. the occurrence of landslides and tsunamis in the near future. Both landslides and tsunamis may represent a danger to offshore infrastructures, associated with present and future hydrocarbons exploitation and other renewable energies (Zhang et al., 2019). Tsunamis may also have an impact on the coastal areas of the nearby regions of northwestern Europe, considering the increasing human pressures of these areas (Imamura et al., 2019). The geological record of the Svalbard continental margin can help us to assess the possible tsunamis induced by future landslides. In fact, the sedimentary record of its continental slope evidences numerous landslides, such as the Storfjorden LS-1 landslide, which forms part of the Storfjorden trough mouth fan, and even other recent landslides located in the interfan area of the Storfjorden and Kveithola trough mouth fans (TMFs) (Pedrosa et al., 2011; Rebesco et al., 2012; Lucchi et al., 2012; Llopart et al., 2015).

The relatively medium-size tsunami potential of Storfjorden LS-1 provides new insights into tsunami wave characteristics and evolution. It will help to better understand possible future tsunami hazard in high latitudes. Geomorphic and geotechnical data have been integrated into the L-ML-HySEA (Landslide Multilayer Hyperbolic Systems and Efficient Algorithms) landslide tsunamigenic model simulating landslide dynamics, tsunami wave generation, propagation and coastal impact.

2 Geological setting

Svalbard is located west of the epicontinental Barents Sea and Norwegian continental margin (Fig. 1a). The archipelago resulted from the opening of the northern Atlantic. The major northwest–southeast fault zone is associated with the Knipovich Ridge and the Hornsund Fault Zone (HFZ), crossing the Spitsbergen Island and spreading western of the Barents (Worsley, 1986; Eiken et al., 1994; Engen et al., 2008; Faleide et al., 2008) (Fig. 1b). The post-rift activity of these fault zones has contributed to deforming the Plio-Pleistocene sedimentary sequence (Faleide et al., 1993, 2008; Fiedler and Faleide, 1996). This fault zone has been reactivated by isostatic loading and unloading rebound periods (Pirli et al., 2013; Newton and Huuse, 2017) that are responsible for the local seismicity in the continental margin and can become the trigger of slope failures (Hampel et al., 2009; L'Heureux et al., 2013; Bellwald et al., 2016). Historical earthquakes with magnitudes up to Mw 5 have been recorded (Auriac et al., 2016) (Fig. 1c).

Figure 1Geographic and geological settings of Storfjorden LS-1. (a) Location map of the study area (blue rectangle). (b) Shaded relief map taken from the International Bathymetric Chart of the Arctic Ocean (IBCAO) version 3.0 (Jakobsson et al., 2012) of the North Atlantic Ocean (Norwegian and Barents seas). The major trough mouth fans (grey polygons) and the major submarine landslides (orange polygons) are located on the map. Compilation from Haflidason et al. (2005), Laberg et al. (2000), Laberg and Vorren (1993), Lindberg et al. (2004), Sejrup et al. (2005), and references therein. KF: Kongsfjorden Fan; IF: Isfjorden Fan; BeF: Bellsund Fan; SF & KvF: Storfjorden Fan and Kveithola Fan; BIF: Bear Island Fan; NSF: North Sea Fan. (c) Shaded relief map of the northwestern Barents Sea (Ottesen et al., 2005) displaying the location of the Hornsund Fault Zone (HFZ) and the historical earthquakes recorded from 1960 to 2018 (source from IRIS catalogue; Incorporated Research Institutions for Seismology). This map also shows the locations of Figs. 2 and 3.

The northwestern Barents continental slope is affected by the Storfjorden and Kveithola TMFs (Fig. 2), created by the high sediment input from the Storfjorden and Kveithola glacial troughs crossing the continental shelf, during the onset of the major Northern Hemisphere glaciations, around 2.6–2.7 Ma (Faleide et al., 1996; Butt et al., 2000; Knies et al., 2009). The Storfjorden TMF seafloor is shaped by the relatively large Storfjorden landslide (LS-1) that extends from the shelf edge to the lower continental slope. In spite of their fresh morphological expression, the seismic stratigraphy indicates that LS-1 (Fig. 3c) is a palaeolandslide above the 0.2 Ma R1 reflector (Rebesco et al., 2012), which is then draped by regional 100 ms ( 8 m, considering 1600 m s−1 sediment velocity) thick sediment units (Lucchi et al., 2012; Llopart et al., 2015), and the sliding mass is a subtabular depositional body, between 3 and 5 ms thick (Pedrosa et al., 2011). The high-resolution seismic stratigraphy in the southern Storfjorden TMF is dominated mainly by an alternation of acoustically stratified and transparent units. Eight seismic stratigraphy units (A to G, from top to bottom) have been identified above the R1 reflector (Fig. 3c). Based on acoustic facies they are as follow: the stratified units A, C, E and G and the transparent units B, D and F. The stratified units are usually continuous with high amplitude and draping the existing topography, while the transparent units present an irregular upper boundary and usually a basal erosive surface that describe individual lenses. Comparable seismic facies have also been found in other TMFs (Laberg and Vorren et al., 1993; Cofaigh et al., 2003).

Figure 2Shaded relief map from merging the regional dataset (colour map, 250 m; Ottesen, 2005) and high-resolution datasets (grey map, 75 m) of the Storfjorden and Kveithola TMFs, where the Storfjorden LS-1 landslide reaches full data coverage.

3 Dataset and methods

3.1 Bathymetric data

For slide modelling, high-resolution multibeam bathymetry datasets from different cruises (SVAIS (The development of an Arctic ice stream-dominated sedimentary system: The southern Svalbard continental margin) onboard BIO (Buque de Investigación Oceanográfica) Hespérides, 2007; EGLACOM (Evolution of a GLacial Arctic COntinental Margin: the southern Svalbard ice stream-dominated sedimentary system) onboard R/V Explora, 2008) have been integrated (Fig. 2). Data processing consisted of cleaning and filtering the navigation data, noise reduction and data editing using the Computer Aided Resource Information System (CARIS) HIPS and SIPS software. Data were gridded at 75 m and partially cover the Storfjorden and Kveithola TMFs ( 15 300 km2). The bathymetric mosaic was completed with lower-resolution bathymetry data at 250 m provided by the Norwegian Hydrographic Service (NHS); these were collected between 1965 and 1985 (Ottesen et al., 2005). Moreover, regional gridded bathymetry data for the Arctic Ocean area (IBCAO version 3.0;, last access: 15 October 2020), interpolated to 2.5 km bin size, were only used for regional figures (Fig. 1b).

Figure 3Storfjorden and Kveithola interfan TMFs. (a) Shaded relief colour map of the southwestern continental slope of the Storfjorden and Kveithola interfan TMFs. The orange line marks the headwall, and the white lines indicate its sidewalls. The pink line corresponds to the multichannel seismic profile (EG-06) acquired during the EGLACOM cruise. (b) Slope gradient values. Artefacts are induced by the slope parallel to the ship tracks. (c) The top corresponds to the Topographic Parametric Sonar (TOPAS) subbottom profile (modified from Llopart et al., 2015) and multichannel seismic profile (EG-06, modified from Rebesco et al., 2014), acquired during the SVAIS cruise. Both profiles are displayed at the same horizontal and vertical scales to show matching of acoustic facies. The bottom shows the line drawing of the multichannel seismic profile displaying the seismic stratigraphy, where the R1 reflector (in blue) is the base of Storfjorden LS-1 and its top is the unit D (Pedrosa et al., 2011; Rebesco et al., 2014; Llopart et al., 2015). TWTT: two-way travel time. (d) Table with the seismic units and subunits, ages, and their lithologies (modified from Llopart et al., 2015).

3.2 Tsunami numerical modelling

The L-ML-HySEA is a mathematical model, which implements a two-phase model to reproduce the interaction between the landslide granular material and the fluid. In the present work, a multilayer non-hydrostatic shallow-water model is considered in order to model the evolution of the ambient water, taking into account dispersive water waves (Fernández-Nieto et al., 2018), and to simulate the kinematics of the Storfjorden LS-1 landslide using the Savage–Hutter model (Eq. 3) (Fernando-Nieto et al., 2008).

The L-ML-HySEA model was validated using laboratory experimental data for landslide-generated tsunamis. A milestone in the validation process of this code consisted in the numerical simulation of the Lituya Bay 1958 mega tsunami with real topobathymetric data obtained from González-Vida et al. (2019). The simulation was also used to generate initial conditions for the Method of Splitting Tsunami (MOST), in order to be initialized for the landslide-generated tsunami scenarios of the National Tsunami Hazard Mitigation Program (NTHMP) mandatory benchmarks in the USA (EDANYA Research Group, HySEA; results for NTHMP's tsunami benchmarking process for landslide-generated tsunamis are available at, last access: 1 January 2015).

L-ML-HySEA needs to incorporate the physical properties of the sediment involved in the landslides. In the Storfjorden (LS-1) case, properties determined by Lucchi et al. (2013) and Llopart et al. (2019) were used. For the purpose of modelling, we have assumed that the landslide took place in a single event. The simulation has been performed by considering a  1.3 critical slope repose angle, since that value has given the best results across the models.

3.2.1 Reconstruction of pre-landslide bathymetry and landslide body geometry

To perform the L-ML-HySEA numerical simulation, it is necessary to reconstruct the pre-bathymetry before the seafloor failure (Macías et al., 2015). For this aim, we used the high-resolution multibeam bathymetry together with seismic profiles published to define the landslide location, its body geometry and its buried thickness (Pedrosa et al., 2011; Lucchi et al., 2012; Rebesco et al., 2012, 2014; Llopart 2015). We assume that the sedimentary infill thickness of 100 ms (Llopart et al., 2015) is roughly similar inside and outside of the landslide and then the present-day bathymetry reproduces the pre-Storfjorden LS-1 100 ms difference.

The pre-Storfjorden LS-1 landslide (Fig. 4a) has been calculated by filling the current headwall and lateral-scarp areas using the cartographic sewing technique on the bathymetry with B splines (Lee et al., 1997) and defining a network of B-spline patches (Eck and Hoppe, 1996). The corresponding control vertices splines were developed using CAD (computer-aided design) software tools through contour lines from the DEM (digital elevation model) and defined by a tolerance rectangle. When creating the spline, the tolerance rectangle is displayed in the form of construction lines. The control vertices of the rectangle, which are shown as circles, influence the spline curves. The spline is tangent to the tolerance rectangle at the start and end points. In this way, the curve ends up adapting to the hypothetical geometry that best fits each patch. Once the splines were developed, patches were densified through the existing DEM and the points were calculated through the splines, generating a new complete DEM without patches. This procedure uses only data points that are not affected by the landslide and assumes convergence of both datasets (boundary conditions) where the slide scars terminate. A second step involves obtaining the volume of the slid sediment body (Fig. 4b) by calculating the difference between the reconstructed pre- and post-landslide bathymetry.

Figure 4The Storfjorden LS-1 landslide geometry. (a) Reconstruction of pre-landslide seafloor morphology. (b) Available displaced sediment thickness in metres. Note the location of the different sliding upper and middle sectors (LS1-U and LS1-M, respectively).

3.2.2 The L-ML-HySEA model equations and discretization

The Multilayer-HySEA model consists of a two-phase model that represents the interaction between a submarine or subaerial landslide (composed by granular material) and the ambient fluid. A multilayer non-hydrostatic shallow-water model is considered for modelling the evolution of the ambient water, and for simulating the kinematics of the submarine/subaerial landslide the Savage–Hutter model is used. As can be found in Fernández-Nieto et al. (2018) the simplest model that can be used for this purpose consists of a single fluid layer (Eq. 1) modelled by a non-hydrostatic shallow-water equation coupled with a Savage–Hutter model (Eq. 2).

(1)S-W systemth+xhu=0,thu+xhu2+12gh2-ghxH-zs=naus-u,(2)S-H systemtzs+xzsus=0,tzsus+xzsus2+12g1-rzs2-g1-rzsxH=-rnaus-u+τP,

where g is the gravity acceleration (g=9.81 m s−2); H(x) is the non-erodible bathymetry measured from a predetermined reference level, zs(x,t) denotes the thickness of the layer of granular material at each point x at time t, h(x,t) is the total water depth, and η(x,t) represents the free surface (measured from the same aforementioned fixed reference level) and is given by η=h+zs-H. u(x,t) and us(x,t) are the averaged horizontal velocity for the water and for the granular material, respectively, and r=ρ1ρ2 is the ratio of densities between the ambient fluid and the granular material. The friction between the fluid and the granular layer is parameterized with the term na(usu). Finally, τP(x,t) represents the friction between the granular slide and the non-erodible bottom surface. The parameterization follows the system proposed in Pouliquen and Forterre (2002).

These two models are coupled through the boundary conditions at their interface. The parameter r represents the ratio of densities between the ambient fluid and the granular material (slide liquefaction parameter).

Usually, it is formulated that

(3) r = ρ f ρ b , ρ b = 1 - φ ρ s + φ ρ f ,

where ρs represents the typical density of the granular material and ρf is the density of the fluid (ρs>ρf), both considered constant, and φ represents the porosity (0φ<1). In this model φ is supposed to be constant in space and time, and, consequently, the ratio r is also constant. This ratio r ranges from 0 to 1 (i.e. 0<r<1) and is a value difficult to estimate even in a uniform material, as it depends on the porosity (and ρf and ρs are also assumed constant) (Fig. 5).

Nevertheless, this system is not accurate enough for modelling landslides where the dispersive effects are relevant, so the model to be considered in this work and described in the next section is the L-ML-HySEA.

Figure 5Schematic figure to describe the multilayer system, with water height (h), unchanged non-erodible bathymetry (H), depth-averaged velocity in the x direction (μα), depth-averaged velocity in the z direction (ωα), sediment thickness (zs), non-hydrostatic pressure at the interface (Pα-1/2), free-surface elevation measured from a fixed reference level (η) and the number of layers (L).


The fluid model

The ambient fluid is modelled by a multilayer non-hydrostatic shallow-water system (Férnandez-Nieto et al., 2018) so that dispersive water waves can be taken into account. The model is obtained by a process of depth averaging of the Euler equations and can be interpreted as a semi-discretization with respect to the vertical coordinate.

The total pressure is decomposed into the sum of hydrostatic and non-hydrostatic components, in order to take into account dispersive effects. In this process, the horizontal and vertical velocities are supposed to have constant vertical profiles. The resulting multilayer model admits an exact energy balance, and when the number of layers increases, the linear dispersion relation of the linear model converges to the same of Airy's theory. Finally, the model proposed in Férnandez-Nieto et al. (2018) can be written in a compact form as

(4) t h + x ( h u ) = 0 , t h u α + x h u α 2 + 1 2 g h 2 - g h x H - z s + u α + 1 2 Γ α + 1 2 - u α - 1 2 Γ α - 1 2 = - h x p α + σ α z p α - τ α , t h w α + x h u α w α + w α + 1 2 Γ α + 1 2 - w α - 1 2 Γ α - 1 2 = - h z p α , x u α - 1 2 + σ α - 1 2 z u α - 1 2 + z w α - 1 2 = 0 ,

for α{1,2,,L}, with L being the number of layers and where the following notation has been used:


where f denotes one of the generic variables of the system (i.e. uα, wα and pα), Δs=1/L and


A schematic picture of model configuration is where the total water height h is decomposed along the vertical axis into L≥1 layers (Fig. 5). The depth-averaged velocities in the x and z directions are written as uα and wα, respectively. The non-hydrostatic pressure at the interface zα+1/2 is denoted by pα+1/2. The free-surface elevation measured from a fixed reference level (for example the still-water level or mean level in the ocean) is written as η and η=h-H+zs, where again H(x) is the unchanged non-erodible bathymetry measured from the same fixed reference level. τα=0, for α>1 and τ1 is given by


where τb stands for a classical Manning-type parameterization for the bottom shear stress and, in this model, is given by


and na(usu1) accounts for the friction between the fluid and the granular layer. The latest two terms are only present at the lowest layer (α=1). Finally, for α=1,,L-1, Γα+1/2 parameterizes the mass transfer across interfaces, and those terms are defined by


Here, we suppose that Γ1/2=ΓL+1/2=0, which means that there is no mass transfer through the seafloor or the water free surface. To close the system, the boundary condition


is imposed at the free surface, and the boundary conditions


are imposed at the bottom. The last two conditions enter into the incompressibility relation for the lowest layer (α=1), given by


It should be noted that the hydrodynamic model described here and the morphodynamic model described in the next subsection are coupled through the unknown zs, which, in the case of the model described here, is present in the equations and in the boundary condition (w0=-tH-zs).

Some dispersive properties of the system (4) were originally studied in (Férnandez-Nieto et al., 2018). Moreover, for a better-detailed study on the dispersion relation (such as “phase velocity”, “group velocity” and “linear shoaling”) the reader is referred to the work of Macías et al. (2020).

Along the derivation of the hydrodynamic model presented here, the rigid-lid assumption for the free surface of the ambient fluid was adopted. Therefore, pressure variations induced by the fluctuation on the free surface of the ambient fluid over the landslide are neglected.

The landslide model

The 1D Savage–Hutter method implemented in the model is given by system (2). The friction law τP (Pouliquen and Forterre, 2002) is given by the expression


where μ is a constant friction coefficient with a fundamental role because it controls the movement of the landslide. Usually, μ is given by the Coulomb friction law as it is the simplest parameterization that can be used in landslide models. However, it is well known that a constant friction coefficient does not allow for models to reproduce the steady uniform flows over rough beds that are observed in the laboratory for a range of inclination angles. In the work of Pouliquen and Forterre (2002), in order to reproduce these flows, the authors introduced an empirical friction coefficient μ that depends on the norm of the mean velocity us on the thickness zs of the granular layer and on the Froude number Fr=usgzs. The friction law is given by




where ds represents the mean size of the grains. β=0.136 and γ=10-3 are empirical parameters. tan (δ1) and tan (δ2) are the characteristic angles of the material, and tan (δ3) is another friction angle related to the behaviour when starting from rest. This law has been widely used in the literature (see for instance Brunet et al., 2017).

It is important to remark that this slide model can also be adapted to simulate subaerial landslides. The presence of the term (1−r), in the definition of the Pouliquen–Forterre friction law, is due to the buoyancy effects, which must be taken into account only in the case that the granular material layer is submerged in the fluid. Otherwise, this term must be replaced by 1 in order to consider subaerial landslides.

In Macías et al. (2021) the reader can find the details about the numerical algorithms used to implement the model.

The discretization of the resulting systems is difficult. For the hydrostatic systems that are expressed as non-conservative hyperbolic systems, the natural extension of the numerical schemes proposed in Escalante et al. (2018, 2019) has been adopted and then solved using a second-order HLL (Harten–Lax–van Leer), positive-preserving, well-balanced, path-conservative finite-volume numerical scheme (see Castro Díaz and Fernandez-Nieto, 2012). Then, the non-hydrostatic pressure corrections at the vertical interfaces required the discretization of an elliptic operator, and that was done using standard second-order central finite differences. This resulted in a linear system that was solved using an iterative scheduled Jacobi method. Finally, the computed non-hydrostatic corrections were used to update the horizontal and vertical momentum equations at each layer, and, at the same time, the frictions were also discretized (see Escalante et al., 2018, 2019). For the discretization of the Coulomb friction term, the procedures presented in Fernández-Nieto et al. (2008) were followed.

The resulting 2D numerical scheme is well balanced for the water at rest stationary solution and is L stable under the normal CFL (Courant–Friedrichs–Lewy) condition. The scheme is also positive preserving; that means that the thickness of the water layer will be always positive or zero but never negative and can be used with emerging topographies.

For dealing with numerical experiments in 2D regions, the computational domain must be decomposed into cells or finite volumes with a simple geometry. Here, a Cartesian-type UTM (Universal Transverse Mercator) was used. The 2D numerical algorithm for the hydrodynamic hyperbolic component of the coupled system is well suited to be parallelized and implemented in GPU (graphics processing unit) architectures, as is shown in Castro et al. (2011). Unfortunately, the standard treatment of the elliptic part of the system is not compatible with the parallelization of the algorithms. However, in Escalante et al. (2018, 2019), a multi-GPU implementation was presented and made possible because of the compactness of the numerical stencil and the massive parallelization of the Jacobi method. Such a multi-GPU implementation of the complete algorithm results in much shorter computational times, and that is the reason why it was used in this work.

4 Results

4.1 The Storfjorden LS-1 landslide geometry

The Storfjorden LS-1 landslide is  60 km in length and covers an area of more than 1300 km2 (Llopart et al., 2015). Three main morphological elements are imaged by the multibeam bathymetry: the headwall, sidewalls and sliding area (Fig. 3a). The headwall displays a well-defined seaward-concave scarp that forms an amphitheatre-like feature about 12 km long and > 50 m in relief. Its slide scar is incised into the shelf edge at 420–480 m water depth. The northwestern sidewall is defined by a striking 25 km long scarp, 35–40 m in relief, with a rectilinear to slightly sinuous pathway. The southeastern sidewall is 35 km long with 25 to 80 m of relief, representing the highest in the mid area ( 1500 m water depth). The width between the sidewalls is variable downslope. The sidewalls are roughly parallel and define a bottleneck shape of 18 km wide, down at  1330 m water depth, which increases to 32 km at  1900 m water depth. The sliding area displays an elongated lobate shape, in plan view, with an irregular seafloor. The seafloor gradients are typically 2 to 3 at  1330 m water depth and < 2 toward the distal ends (Fig. 3b). The chaotic landslide deposits in the upper slope (800 m depth) are shown in Fig. 3c (Llopart et al., 2015).

4.2 Submarine landslide and tsunami numerical simulations

The numerical simulation consists of several successive steps aimed at reconstructing (i) the smooth pre-landslide upper slope and landslide body geometry following the methods described in Sect. 3.2.1 (Fig. 4), (ii) the landslide dynamic (Fig. 6), (iii) the tsunami wave generation (Fig. 7), and (iv) the tsunami wave propagation and its impacts on the coast (Figs. 8 and 9).

4.2.1 Modelling the landslide dynamic

Once the smooth pre-landslide upper slope had been calculated, following the methods described in Sect. 3.2.1 (Fig. 4), the landslide body geometry was determined. The numerical landslide rupture simulation begins with the slope failure of Storfjorden LS-1, which assumes that it fails at once and moves downslope by gravitational forces (Macías et al., 2015). Conventional studies about submarine slides show the difficulty in assessing whether their occurrence represents unique events (Casas et al., 2016; Chiocci and Casalbore, 2017; Vázquez et al., 2022). The morphological results would also support that assumption for Storfjorden LS-1, due to the lack of retrogressive structures, which would point to a decrease in the tsunamigenic potential (Harbitz et al., 2006). A few local and small-scale slope failures seem to occur on the southeastern flank in the mid slope ( 1600 m deep), but they would not be significant in the modelling of the main landslide event, due to its low potential to transfer deformation to the water column.

In addition, by comparison with similar deposits from other continental slopes (Iglesias et al., 2012; Casas et al., 2013; Casas et al., 2013; Vanneste et al., 2006; Winkelmann et al., 2008), Storfjorden LS-1 shows a well-defined arcuate slide scar, mostly rectilinear side walls and a cutting basal shear surface. Its moving mass defines a subtabular body with chaotic deposits and without apparent internal discontinuities that sharply interrupt the lateral continuity of the surrounding deposits (Pedrosa et al., 2011; Lucchi et al., 2012). All these characteristics put together tentatively point to a single sliding process, not being demonstrated as multiple failures in the previous literature (Pedrosa et al., 2011; Lucchi et al., 2012; Rebesco et al., 2012, 2013). The difference between pre- and post-landslide bathymetries is 40 km3, a volume between 33 km3 proposed by Pedrosa et al. (2011) and 46 km3 by Llopart et al. (2015). The numerical landslide rupture simulation shows that the moving mass was comprised of two domains with different behaviour, based on the velocity pattern (Fig. 6a and b). At 1 min after the slide, an upper-slope domain (LS1-U,  500 m water depth,  3 slope gradients) is related to the moving sediment nearest to the slide scar and moves faster (vu=5 m s−1) than the mid-slope domain (LS1-M,  1000 m water depth, 2 slope gradients, vm=1 m s−1) (Fig. 6a, step 1; Fig. 7, step 1; Supplement video 1, EDANYA Research Group, 2022a). At 4 min, the velocity increases in the upper-slope domain as vu= 30 m s−1 and in the mid-slope domain as vm=22 m s−1 (Fig. 6, step 2; Fig. 7, step 2; Supplement video 1). At 20 min, the sliding velocity maximum in the frontal area is vf= 20 m s−1 and decreases gradually towards the sides of the slide with values around vs= 15 m s−1 (Fig. 6c, step 10; Fig. 7, step 10; Supplement video 1). At 25 min, the sediment sliding became homogeneous in the distal area ( 1500 m water depth). The velocity decreases gradually downslope, where the frontal area of sediment sliding velocities reach around vf= 10 m s−1 and in distal area vs= 2–5 m s−1.

The landslide characteristics modelled by L-ML-HySEA determine an average velocity (va) of 25 m s−1, a terminal time (tt) of  40 min and a characteristic distance (dc) of at least  60 km.

Figure 6Velocity pattern of the landslide in different frames and the displaced sediment thickness. (a) A 1 min frame, where the upper-slope domain (LS1-U,  500 m water depth, 3 slope gradient) with an average velocity vu= 5 m s−1 contrasts with the mid-slope domain (LS1-M,  1000 m water depth; 2 slope gradients) with an average velocity vm= 1 m s−1. (b) A 4 min frame, with a progressive velocity increase in both areas (vu= 30 m s−1 and vm= 22 m s−1, respectively). (c) A 20 min frame, where velocity is maximum in the frontal area as vf= 20 m s−1, coinciding with maximum sediment thickness displacement (122–108 m) and decrease towards the sides as vs= 15 m s−1. (d) A 25 min frame, where both moving masses are joined with central velocity vf= 10 m s−1 and lateral velocity vs= 2–5 m s−1 in the distal area ( 1500 m water depth).

Table 1Storfjorden LS-1 geometry and mechanical characteristics. The inputs for the model are related to the landslide geometry: the total volume (V), the longitude and latitude of the submarine landslide relative to the LS1-U slope area (X1, Y1) and the LS1-M slope area (X2, Y2), the initial depth (h) before the slope failure, the length (l) (long axis) and width (w) (small axis), the maximum thickness in the proximal and distal areas (Tp and Td, respectively), the mean azimuth direction of the landslide (Az), and the mean slope gradient (O). The outputs of the model are related to the landslide dynamic: velocities to the landslide in the upper slope (vu) and the landslide located in the mid slope (vm), sliding sediment in the frontal area (vf), sliding sediment in the sides (vd), total velocity (Vt), terminal duration (tt), and characteristic distance (dc). The outputs related to initial tsunamis: wave crest (Ci) tsunami velocity at 1900 m depth (vts), tsunami wave velocity during refraction (vr), tsunami wave velocity toward the north (vn) and tsunami wave velocity toward the southeast (ve). LS1-U: Storfjorden landslide upper slope; LS1-M: Storfjorden landslide mid slope.

Download Print Version | Download XLSX

4.2.2 Tsunami wave generation

Free-water-surface changes at any point along the time are determined according to seafloor deformation (Fig. 7 and Supplement video 1). In the initial stage (4 to 5 min), the tsunami wave has two northwest–southeast-trending dipoles, LS1-U, smaller (25 km long) and more striking, and LS1-M, larger (35 km long) and smoother. They have been created by the water mass infilling the empty spaces produced by the sudden evacuations and uplifting of fast-downslope-moving mass. Both wave dipoles have the troughs in shallower waters ( 600 to 800 m) than their respective crests ( 700 to 1000 m water depth) (Fig. 7, steps 1 and 8). The synthetic marigram on the upper slope (station 1) of Storfjorden LS-1 (Fig. 9a, b and c) highlights the initial wave generation with the crest and trough well defined, registering a crest amplitude value of 0.4 m (above LS1-U) and a trough amplitude value of up to 1.1 m at 3 min, followed by a crest amplitude value of 0.3 m (Fig. 8, step 1; Fig. 9b; Supplement video 2, EDANYA Research Group, 2022b).

After 3 min, the two initial dipoles evolve into a single northwest–southeast-trending dipole (crest amplitude value of 0.3 m), whose trough (0.5 m) is also in shallower waters than the crest (Fig. 7, step 3; Fig. 8, step 1). Wave rebound occurs (Fig. 7, steps 4 to 9) when a maximum amplitude of  2.7 m is registered over the distal area of the landslide (Fig. 9c, station 2; Supplement video 2). At 7 to 9 min, a new crest wave amplitude value of 0.5 m appears parallel to the single dipole (Fig. 7, steps 4 to 5; Supplement video 1) in shallower waters. It enlarges with time up to 0.7 m, whereas the trough largely keeps its dimensions or relief. At the 16 min mark, the tsunami wave reaches the highest amplitude (at 1780 m water depth), with a trough amplitude value of 2 m and a crest amplitude value of 0.7 m (Fig. 9, station 2).

Figure 7Main composition in 12 consecutive frames. (a) Evolution of the landslide, where two sliding sediment masses in the firsts 25 min can be distinguished, one in the upper slope (LS1-U) and the other one in the mid slope (LS1-M), which are merged into a unique one after 25 min. The colour scale corresponds to available sediment thickness to be displaced in metres. (b) Dipole wave evolution of the outgoing tsunami with the wave height values.

At 25 min, the tsunami wave evolves into a larger dipole above the landslide, opposite to the first ones, with a trough amplitude value of 0.5 m and a crest amplitude value of 0.3 m (at 1200 m water depth). This dipole gets smaller with time, with a crest amplitude value of 0.5 m in shallower waters (900 m water depth) covering large areas with time (Fig. 7, step 10).

4.2.3 Tsunami wave propagation and coastal impact

The tsunami wave dynamics are illustrated by the maps of the wave height across time (Fig. 8 and Supplement videos 1 and 2). Synthetic marigrams have been included at key locations in order to highlight the wave propagation and coastal impact in the northwestern flank of Spitsbergen Bank (Fig. 9d, station 3), the onshore area of the Kveithola glacial trough (Fig. 9e, station 4), Spitsbergen Bank (Fig. 9f, station 5), the onshore northern boundary of the mid shelf of Storfjorden glacial trough (Fig. 9g, station 7) and the onshore southwestern Spitsbergen coast (Fig. 9h and i, stations 7 to 8).

The initial tsunami wave starts propagating from the landslide area with a trough wave moving northeast towards the coast and a crest wave moving southwest (Fig. 8, step 3; Supplement videos 1 and 2). The tsunami waves propagate elliptically, with the crest and trough elongated in the northwest–southeast direction (Fig. 7, step 3) and a total velocity (vts) of  136 m s−1. During the tsunami propagation from depth ( 1900 m) towards shallow water ( 250 m) (Fig. 8 and Supplement video 2), the refraction phenomenon occurs during shoaling. The synthetic marigram records the refraction clearly due to an irregular variation in the general amplitude pattern (Fig. 9d, station 3). Refraction effects in a wave front show a decrease in velocity down to vr= 81 m s−1 and an increase in amplitude. The values change from a trough amplitude of 0.25 m to a crest amplitude value of 0.18 m (Fig. 9e, station 4), and the values change from a trough amplitude of 4.2 m to a crest amplitude of 4.3 m (Fig. 9f, station 5). Furthermore, it can be observed that the tsunami propagation front displays a crescent shape.

Figure 8Eight steps displaying the tsunami waves generated at different times. The dotted red line indicates the shelf break. The first step shows the generation of the tsunami wave and its spread. At 23 min the outgoing wave suffers refraction. At 50 min the tsunami wave arrives at the northern coast of Spitsbergen and Bear Island, and at the 60 min, it hits the westernmost coast of Spitsbergen.

The tsunami arrival times are observed based on the propagation direction toward the coast. The impact of tsunami waves affects Sørkappøya at 50 min (southern Spitsbergen), with a trough amplitude value of 0.3 m, increasing to a crest amplitude value of 0.2 m at 75 min (Fig. 8, step 5; Fig. 9h). At the same time (50 min), the impact occurs at Kapp Dunér (northwestern Bear Island) with tsunami waves having a crest amplitude value of 0.3 m at 50 min, which increases to 0.5 m at 53 min (Fig. 8, step 5; Supplement video 2). After these two first coastal impacts, the tsunami affects different parts of both islands at different times. Southwestern Bear Island is reached by the tsunami waves at 60 min with a maximum crest amplitude value of 0.5 m (Fig. 8, step 6), followed by a trough amplitude value of 0.5 m at 65 min. Likewise, the north of Bear Island is affected by a trough amplitude value of 0.5 m (Fig. 8, step 6) and a crest amplitude value of 0.5 m at 80 min (Fig. 8, step 7). In Svalbard, Stormbukta Bay (southwestern Spitsbergen Island) is the next impacted coastline (Fig. 9a). There, the tsunami waves show specific velocity of vs= 13 m s−1 (at 18 m of depth), with trough amplitude values of 0.3 m at 63 min that increase crest amplitude values up to 0.32 m (85 to 95 min) (Fig. 8, steps 7 to 8; Fig. 9i). Finally, tsunami wave series propagating toward the coast occur until 2 h after the Storfjorden LS-1 landslide triggering (Supplement video 2).

Figure 9Synthetic marigrams. (a) The location map of the synthetic marigrams above the reference stations, numbered in panel (a) (station 1 to 8). (b–i) Note the different arrival times, periods, wave heights and polarities at the different stations.

5 Discussion and conclusions

The L-ML-HySEA landslide tsunamigenic model provides a fast and consistent method for simulating landslide dynamics, tsunami wave generation, propagation and coastal impact. The results reveal several fundamental insights regarding the assessment of the main factors that control the characteristics and evolution of the Storfjorden LS-1 tsunamigenic landslide, as well as the related coastal hazard.

The submarine landslide geohazard studies are not easy to conduct due to the difficult access that the marine environment imposes, which makes it hard to obtain accurate analysis and monitoring. Our case study suggests that numerical modelling is of great help in understanding the dynamics of submarine landslide geohazard and their tsunamigenic potential. The Storfjorden upper continental slope presents critical conditions that need to be taken into account and warrants carrying out studies to assess slope stability. Several factors support this assertion: (i) the overpressure ratios measured in the subsurface sediments (Lucchi et al., 2013; Llopart et al., 2019) and (ii) the seismicity related to the active Hornsund Fault Zone (Hampel et al., 2009; Auriac et al., 2016; Pirli et al., 2013). In addition, (iii) the recent environmental stress represented by the factors mentioned above may intensify processes such as gas hydrate dissociation and fluid flow migration (León et al., 2021). Lastly, (iv) unloading rebound seismicity during ice retreat/melting (Berndt et al., 2009) also may contribute to triggering new submarine landslides.

5.1 Landslide dynamics and wave generation

Key landslide parameters for the generation of tsunamis commonly include volume, velocity and initial acceleration of the sliding mass (Harbitz et al., 2006; Urlaub et al., 2013; Løvholt et al., 2015; Macías et al., 2015; Urgeles et al., 2019) (Table 1). In our case study, the relationship between the above-mentioned parameters is also fundamental for the formation of a tsunami. The L-ML-HySEA model indicates that the 40 km3 of available displaced sediment volume moved by the Storfjorden LS-1 landslide is enough to trigger a tsunami. The morphosedimentary characteristics of the Storfjorden LS-1 landslide suggest that mass failure deposits could occur as a single event, as opposed to several, implying a better energy transfer to the water column (Vázquez et al., 2022, and references therein); therefore, velocity and initial acceleration also would be key to the formation of the tsunami. In this sense, for the phase velocity to be highly effective at the depths of Storfjorden LS-1 (H= 420 to 1900 m) during the tsunami wave onset, its value should be vta= 136 m s−1 (Tinti et al., 2000; Fryer et al., 2004). The relatively average velocity (roughly 25 m s−1, Table 1) obtained for our tsunami indicates that it was out of phase, and, therefore, it would not have been effective enough to create high-amplitude tsunami waves (Huggel et al., 2005; Evans et al., 2009; Pudasaini, 2014; Dietrich and Krautblatter, 2019).

Our results indicate that the characteristics of the tsunami wave are influenced by landslide dynamics at two stages of the downslope moving mass: initial (1 to 4 min) and late (20 to 25 min), with different velocity values (Fig. 6). The two initial wave dipoles, generated when seafloor failure occurs, are the consequence of two large seafloor depressions, one at the slide scar (i.e. main evacuation area) and the other one located between the two mass moving domains (Fig. 7 and Supplement video 1). At 25 min, a new dipole is formed, opposite to the previous ones. This new dipole is created when the faster sliding mass reaches the slower one and both masses merge producing a significant impact in the available displaced sediment thickness (100 m) of the distal moving mass. The increase in thickness would contribute to increasing the pressure in the water column, causing the uplift in the water surface and the enhancement of the tsunami waves (Ramadan et al., 2018; Ercilla et al., 2021). Thus, our study demonstrates that a proper understanding of landslide dynamics at their initial stages (or first motion) and of their deformation during the run-out is a crucial requirement for understanding the characteristics of the initial tsunami waves and the effects that those characteristics have on their evolution. In addition, our study also suggests that identifying the initial tsunami wave forms could reveal the tsunami sources, e.g. landslides (generating single or multiple trough and crest pairs) versus faults (generating a single or crest wave) (e.g. Macías et al., 2015; Ercilla et al., 2021; Estrada et al., 2021; Bécel et al., 2017).

5.2 Seafloor morphology

It is widely known that how a tsunami wave propagates is highly dependent on the morphology of the seafloor (Urlaub et al., 2013; Estrada et al., 2021). The model shows that tsunami waves propagate elliptically with respect to the northeast–southwest elongated seafloor shape of the landslide relief. In the Storfjorden LS-1 tsunami, the average velocity waves (vta= 136 m s−1) travel northeastern and are focused between the continental shelf of Svalbard and of Bear Island, which helps to confine it and forces its direction of propagation (Fig. 8 and Supplement video 2). The continental shelf morphology determines tsunami shoaling, with the shallowest water depths located at Spitsbergen Bank (80 m water depth). The shoaling by the bank produces the refraction phenomenon and the amplification of the tsunami. As the tsunami propagates across the ocean, waves can undergo refraction, which is caused by segments of the wave moving at different speeds as the water depth along the wave front varies (Berkhoff, 1972). This effect produces variations in the amplitude of the tsunami waves and in their arrival times: 15 min later at southwestern Spitsbergen and 11 min later at northwestern Bear Island. When encountering an obstacle, the tsunami waves discharge their energy with great force, as in the case of Sørkappøya and Sørkapp at 75 min, decreasing the amplitude and slowing down the arrival time of the tsunami waves in the corresponding bay to 100 min (Supplement video 2). The tsunami wave arrival is recorded at 80 to 95 min at the eastern Hornsund fjord. On the other hand, the crescent shapes of the tsunami front seem to be conditioned by the Storfjorden and Kveithola glacial troughs separated by Spitsbergen Bank. The elongated negative reliefs of the glacial troughs would cause the funnelling of the tsunami seawater with relative higher specific velocities of propagation (vs= 56 m s−1 in the Storfjorden glacial trough at 320 m water depth). Therefore, the numerical simulations are a useful tool to assess tsunami hazard in places where local seafloor topography could advance or delay the tsunami waves and therefore the coastal impact.

The Storfjorden LS-1 modelling has demonstrated that the shelf seafloor morphology is a decisive factor: it influences the propagation velocity of the tsunami waves, the variations in wave amplitude (shoaling effect) and the impact of coastal arrivals (Iglesias et al., 2012; Estrada et al., 2021; Salaree and Okal, 2020).

5.3 Coastal location

The initial tsunami waves start propagating as negative and positive disturbance dipoles. The trough is always located towards the upper part of the margin, which determines that, in general, the first arrival to the coast corresponds to a sea level drop, hence decreasing the coastal impact. This factor, together with coastal location and orientation (i.e. angle of the waves with the coastline), conditions the polarities of the first arrival wave. In our study, the first arrival wave impacted with trough values of 0.1 to 0.5 m against the coast of southwestern Spitsbergen and crest values of 0.1 to 0.5 m against Bear Island (Fig. 8 and Supplement video 2).

5.4 Comparison with other landslide tsunamis: hazard assessment for the northern glaciated margins

Landslide parameters (age, area, volume, seafloor gradients, location and velocity) and related tsunami parameters (wave amplitude and velocity) of different landslide-inducing tsunamis have been compared with those defining Storfjorden LS-1 and its tsunami, in Table 2.

Table 2Submarine landslides compared with Storfjorden LS-1. In the trigger column, the letters represent (A) weak-layer sedimentary architecture, (b) earthquake, (c) gas hydrate dissociation and (d) other additional factors. Deposit type: mass wasting (MW), large rafted blocks (LRBs), turbidity currents (TCs), debris flows (DFs) and slump blocks (SBs). References: Storfjorden LS-1 (Pedrosa et a., 2011; Lucchi et al., 2012; Rebesco et al., 2014; Llopart et al., 2015), Hinlopen (Winkelmann et al., 2008; Vanneste et al., 2006), Bjørnøyrenna (Laberg et al., 1999), Gebra Valley (García et al., 2009; Casas et al., 2013), Grand Banks (Piper et al., 1999; Fine et al., 2005), Storegga (Haflidason et al., 2004; Bondevik et al., 2012), Kongsfjorden TMF (Berndt et al., 2009), Trænadjupet (Laberg et al., 2002), Nyk (Lindberg et al., 2004), Skagway (Thomson et al., 2001; Synolakis et al., 2002), Lituya Bay (González-Vida et al., 2019) and BIG'95 (Iglesias et al., 2012).

Download Print Version | Download XLSX

This comparison highlights that not only large but also relatively small–medium-sized landslides could have triggered tsunamis in the past. Despite Storfjorden LS-1 having a smaller area and volume with respect to the larger landslides of the Kongsfjorden TMF, the amplitude value of their respective tsunami waves is roughly similar. The effects of global warming over the landslide-triggering factors (e.g. isostatic rebound seismicity and gas hydrate destabilization by the rise in temperature of the ocean water) is not likely to provoke the occurrence of such large landslides as those formed during glacial maxima and the transition from glacial to interglacial periods (Lee, 2009, and references therein). However, the present trend of global warning should over time increase the probability of slope instability, especially on those glaciated margins that have not yet failed after the last glacial to interglacial transition, for instance, the Bear Island, Kongsfjorden and Storfjorden TMFs (Berndt et al., 2009). Therefore, the results presented here should encourage us to continue working in the prediction of tsunamigenic landslide hazards and their coastal impact, mainly in the northern glaciated margins.

The water depth location of the landslide scar seems to influence the coastal impact of tsunami waves (Table 2). In the study area, they present low-amplitude values at the coastal area and their arrival times are longer (50 to 80 min) than the tsunamis modelled in the nearby coast of western Spitsbergen (Bernt et al., 2009), where the slope failure is at shallower water depths (200 m) and the distance to the nearby coast is shorter ( 90 km). Landslides triggered in shallower water result in more localized waves, and the elongated landslide velocity profile delays the appearance of the first positive landward-propagating wave, hence reducing the chances of constructive interference along the coast (Harbitz et al., 2006). This suggests that tsunami modelling based on past landslides should pay more attention to those sectors of the northern glaciated margins with narrower continental shelves and submarine landslide head scarps that are near the shoreline. The landslides located in the middle and northern parts of Norway are a good example of that (L'Heureux et al., 2013; Beaten et al., 2013). Also, much smaller collapses, either submarine or subaerial, also pose a significant local threat. In Norway, several rockslide tsunamis occurred in the 20th century, the most devastating in Storfjorden in 1934 (Blikra et al., 2005; Böhme et al., 2015). Although this is not the scope of the paper, they could be more frequent with an estimated rate of 1 event per 1000 years (Blikra et al., 2005).

In summary, our findings demonstrate that tsunami modelling based on past landslides using the L-ML-HySEA landslide tsunami model will be useful to provide new perspectives on tsunami hazard assessment in polar margins, where global climatic change and its related ocean warming may contribute to the activation of landslides. Landslide tsunami models will allow us to identify the areas with maximum and faster coastal impact and the effect of the local bathymetry on tsunami direction of propagation, shoaling, amplification and diffraction. This knowledge is very important for the design of early-warning strategies, as it will contribute to assess the key factors that are useful as emergency planning tools.

Code and data availability

The SVAIS data are hosted by the Oceanography Cruise and Data Catalogue of CSIC-UTM/Marine Technology Unit (CSIC, 2018).

Video supplement

Supplement video 1 (EDANYA Research Group, 2022a), showing the Storfjorden landslide tsunami, is available at Supplement video 2, showing Storfjorden landslide tsunami propagation, is available at (EDANYA Research Group, 2022b).

Author contributions

MTPG, JMGV and GE conceived the idea for the study and wrote the paper with methodological contributions from JMGV in the tsunami mathematical model and SO in the seafloor deformation and the tsunami potential. The geodynamic framework was done by MTPG, JMGV and GE. Seismic data compilation and the figures were done by MTPG. In the same manner DC contributed to the discussion. The propagation tsunami model was done by JMGV, SO and MJC, as well as associated videos. All authors provided guidance on the analyses and commented on the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Tsunamis: from source processes to coastal hazard and warning”. It is not associated with a conference.


We thank Dimitris Sakellariou and Daniele Casalbore for their helpful comments and discussion. This research was supported by Spanish IPY projects SVAIS (no. POL2006-07390/CGL) and IPY-NICE STREAMS (no. CTM2009-06370-E/ANT; Neogene ice streams and sedimentary processes on high-latitude continental margins, incorporated into the International Polar Year as activity no. 367), as well as the research group RNM 148 and projects PAPEL (no. B-RNM-301-UGR18) and AGORA (no. P18-RT-3275) (Junta de Andalucía, FEDER and University of Granada) and GOLETA (no. PID2019-108880RJ-I00/AEI/10.13039/501100011033). We would like to thank the EDANYA (Differential Equations, Numerical Analysis and Applications) research group from the Universidad de Málaga and also Javier Valencia for the reconstruction of the pre-landslide bathymetry. We also thank all scientists and crew who participated in seagoing activities to obtain geophysical data in the framework of the SVAIS cruise. We are grateful to Dag Ottesen (Norges Geologiske Undersøkelse, NGU) for providing the Norwegian Hydrographic Service bathymetry. The authors wish to acknowledge the cooperation of Captain Pedro Luis de la Puente García-Ganges (BIO Hespérides) and of the technical staff at the UTM (Unidad Tecnología Marina; Consejo Superior de Investigaciones Científicas, CSIC, Barcelona). This work represents a contribution to the CSIC Thematic Interdisciplinary Platform PTI POLARCSIC.

Financial support

This research has been supported by the Consejería de Universidad, Investigación e Innovación, Junta de Andalucía (grant nos. B-RNM-301-UGR18, P18-RT-3275 and RNM 148; University of Granada, FEDER) and the Agencia Estatal de Investigacíon (grant nos. PID2019-108880RJ-I00/AEI/10.13039/501100011033, POL2006-07390/CGL and CTM2009-06370-E/ANT).

We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Review statement

This paper was edited by Hélène Hébert and reviewed by Dimitris Sakellariou and Daniele Casalbore.


Auriac, A., Whitehouse, P. L., Bentley, M. J., Patton, H., Lloyd, J. M., and Hubbard, A.: Glacial isostatic adjustment associated with the Barents Sea ice sheet: a modelling inter-comparison, Quaterny Sci. Rev., 147, 122–135,, 2016. 

Bécel, A., Shillington, D. J., Delescluse, M., Nedimović, M. R., Abers, G. A., Saffer, D. M., and Kuehn, H.: Tsunamigenic structures in a creeping section of the Alaska subduction zone, Nat. Geosci., 10, 609–613,, 2017.  

Bellwald, B., Hjelstuen, B. O., Sejrup, H. P., and Haflidason, H.: Postglacial mass movements and depositional environments in a high-latitude fjord system–Hardangerfjorden, Western Norway, Mar. Geol., 379, 157–175,, 2016. 

Berkhoff, J. C. W.: Computation of combined refraction – diffraction, Coastal Engineering Proceedings, 1, 471–490,, 1972. 

Berndt, C., Brune, S., Nisbet, E., Zschau, J., and Sobolev, S. V.: Tsunami modelling of a submarine landslide in the Fram Strait, Geochem. Geophys., 10, Q04009,, 2009. 

Blikra, L.H., Longva, O., Harbitz, C., and Løvholt, F.: Quantification of rock-avalanche and tsunami hazard in Storfjorden, western Norway, in: Landslides and Avalanches, ICFL Norway, edited by: Senneset, K., Flaate, K., and Larsen, J. O., Taylor and Francis, London, IK, 57–64, ISBN 0415386780, 2005. 

Böhme, M., Oppikofer, T., Longva, O., Jaboyedoff, M., Hermanns, R. L., and Derron, M. H: Analyses of past and present rock slope instabilities in a fjord valley: Implications for hazard estimations, Geomorphology, 248, 464–474,, 2015. 

Bondevik, S., Mangerud, J., Dawson, S., Dawson, and A., and Lohne, Ø.: Record-breaking height for 8000-year-old tsunami in the North Atlantic, EOS, 84, 289–297,, 2003. 

Bondevik, S., Stormo, S. K., and Skjerdal, G.: Green mosses date the Storegga tsunami to the chilliest decades of the 8.2 ka cold event, Quaterny Sci. Rev., 45, 1–6,, 2012. 

Brunet, M., Moretti, L., Le Friant, A., Mangeney, A., Fernández Nieto, E. D., and Bouchut, F.: Numerical simulation of the 30–45 ka debris avalanche flow of Montagne Pelée volcano, Martinique: from volcano flank collapse to submarine emplacement, Nat. Hazards, 87, 1189–1222,, 2017. 

Bugge, T., Befring, S., Belderson, R. H., Eidvin, T., Jansen, E., Kenyon, N., Holtedahl, H., and Sejrup, H. P.: A giant three-stage submarine slide off Norway, Geo-Mar. Lett., 7, 191–198, 1987. 

Butt, F. A., Elverhøi, A., Solheim, A., and Forsberg, C. F.: Deciphering late Cenozoic development of the western Svalbard Margin from ODP Site 986 results, Mar. Geol., 169, 373–390,, 2000. 

Canals, M., Lastras, G., Urgeles, R., Casamor, J. L., Mienert, J., Cattaneo, A., De Batist, M., Haflidason, H., Imbo, Y., Laberg, J. S., Locat, J., Long, D., Longva, O., Masson, D. G., Sultan, N., Trincardi, F., and Bryn, P.: Slope failure dynamics and impacts from seafloor and shallow sub-seafloor geophysical data: case studies from the COSTA project, Mar. Geol., 213, 9–72,, 2004. 

Casas, D., Ercilla, G., García, M., Yenes, M., and Estrada, F.: Post-rift sedimentary evolution of the Gebra Debris Valley. A submarine slope failure system in the Central Bransfield Basin (Antarctica), Mar. Geol., 340, 16–29,, 2013. 

Casas, D., Chiocci, F., Casalbore, D., Ercilla, G., and De Urbina, J. O.: Magnitude-frequency distribution of submarine landslides in the Gioia Basin (southern Tyrrhenian Sea), Geo-Mar. Lett., 36, 405–414, 2016. 

Castro, M. J., Ortega, S., De la Asuncion, M., Mantas, J. M., and Gallardo, J. M.: GPU computing for shallow water flow simulation based on finite volume schemes, C. R. Mec., 339, 165–184,, 2011. 

Castro Díaz, M. J., and Fernández-Nieto, E.: A class of computationally fast first order finite volume solvers: PVM methods, SIAM, J. Sci. Comput., 34, A2173–A2196,, 2012. 

Chiocci, F. L. and Casalbore, D.: Unexpected fast rate of morphological evolution of geologically-active continental margins during Quaternary: Examples from selected areas in the Italian seas, Mar. Petrol. Geol., 82, 154–162, 2017. 

Chiocci, F. L. and Ridente, D.: Regional-scale seafloor mapping and geohazard assessment. The experience from the Italian project MaGIC (Marine Geohazards along the Italian Coasts), Mar. Geophys. Res., 32, 13–23,, 2011. 

Cofaigh, C. Ó., Taylor, J., Dowdeswell, J. A., and Pudsey, C. J.: Palaeo-ice streams, trough mouth fans and high‐latitude continental slope sedimentation, Boreas, 32, 37–55,, 2003. 

CSIC: Oceanography Cruise and Data Catalogue of CSIC-UTM/Marine Technology Unit, CSIC [data set],, 2018. 

Dietrich, A. and Krautblatter, M.: Deciphering controls for debris-flow erosion derived from a LiDAR-recorded extreme event and a calibrated numerical model, Earth Surf. Proc. Land., 44, 1346–1361,, 2019. 

Dowdeswell, J. A., Ottesen, D., Evans, J., Cofaigh, C. Ó., and Anderson, J. B.: Submarine glacial landforms and rates of ice-stream collapse, Geology, 36, 819–822,, 2008. 

Eck, M. and Hoppe, H.: Automatic reconstruction of B-spline surfaces of arbitrary topological type, in: Proceedings of the 23rd annual conference on Computer Graphics and Interactive Techniques, 325–334,, 1996. 

EDANYA Research Group: Video 1 – Storfjorden landslide tsunami, TIB AV-Portal [video],, 2022a. 

EDANYA Research Group: Video 2 – Storfjorden landslide tsunami propagation, TIB AV-Portal [video],, 2022b. 

Eiken, O.: Seismic Atlas of Western Svalbard, Norsk Polarinstitutts Meddelelser, 130, 73, ISBN 82-7666-067-3, 1994. 

Elverhøi, A., De Blasio, F. V., Butt, F. A., Issler, D., Harbitz, C., Engvik, L., and Marr, J.: Submarine mass-wasting on glacially-influenced continental slopes: processes and dynamics. Geological Society, London, Special Publications, 203, 73–87, 2002. 

Engen, Ø., Faleide, J. I., and Dyreng, T. K.: Opening of the Fram Strait gateway: A review of plate tectonic constraints, Tectonophysics, 450, 51–69,, 2008. 

Ercilla, G., Casas, D., Alonso, B., Casalbore, D., Galindo-Zaldívar, J., García-Gil, S., Martorelli, E., Vázquez, J. T., Azpiroz-Zabala, M., DoCouto, D., Estrada, F., Fernández-Puga, M. C., González-Castillo, L., González-Vida, J.M., Idárraga-García, J., Juan, C., Macías, J., Madarieta-Txurruka, A., Nespereira, J., Palomino, D., Sánchez-Guillamón, O., Tendero-Salmerón, V., Teixeira, M., Valencia, J., and Yenes, M.: Offshore Geological Hazards: Charting the Course of Progress and Future Directions, Oceans, 2, 393–429,, 2021. 

Ercilla, G., Casas, D., Alonso, B., Casalbore, D., Estrada, F., Idárraga-García, J., López-González, N., Pedrosa, M., Teixeira, M., Sánchez-Guillamón, O., Azpiroz-Zabala, M., Bárcenas, P., Chiocci, F. L., García, M., Galindo-Zaldívar, J., Geyer, A., Gómez-Ballesteros, M., Juan, C., Martorelli, E., and Yenes, M.: Deep Sea Sedimentation, edited by: Jack, J. and Shroder, F., in: Treatise on Geomorphology (Second Edition), Academic Press, 960–988,, ISBN 9780128182352, 2022. 

Escalante, C., Morales de Luna, T., and Castro, M. J.: Non-hydrostatic pressure shallow flows: GPU implementation using finite volume and finite difference scheme, Appl. Math. Comput., 338, 631–659,, 2018. 

Escalante, C., Dumbser, M., and Castro, M. J.: An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous Galerkin schemes, J. Comput. Phys., 394, 385–416,, 2019. 

Estrada, F., González-Vida, J. M., Peláez, J. A., Galindo-Zaldívar, J., Ortega, S., Macías, J., and Ercilla, G.: Tsunami generation potential of a strike-slip fault tip in the westernmost Mediterranean, Sci. Rep.-UK, 11, 1–9,, 2021. 

Evans, S. G., Bishop, N. F., Smoll, L. F., Murillo, P. V., Delaney, K. B., and Oliver-Smith, A.: A re-examination of the mechanism and human impact of catastrophic mass flows originating on Nevado Huascarán, Cordillera Blanca, Peru in 1962 and 1970, Eng. Geol., 108, 96–118,, 2009. 

Faleide, J. I., Solheim, A., Fiedler, A., Hjelstuen, B. O., Andersen, E. S., and Vaneste, K.: Late Cenozoic evolution of the western Barents Sea-Svaldbard continental margin, Global Planet. Change, 12, 53–74,, 1996. 

Faleide, J. I., Vågnes, E., and Gudlaugsson, S. T.: Late Mesozoic–Cenozoic evolution of the southwestern Barents Sea in a regional rift-shear tectonic setting, Mat. Pet. Geol., 10, 186–214,, 1993. 

Faleide, J. I., Tsikalas, F., Breivik, A. J., Mjelde, R., Ritzmann, O., Wilson, J., and Eldholm, O.: Structure and evolution of the continental margin off Norway and the Barents Sea, J. Int. Geosci., 31, 82–91,, 2008. 

Fernández-Nieto, E. D., Bouchut, F., Bresch, D. Castro., and M. J., Mangeney, A.: A new Savage–Hutter type model for submarine avalanches and generated tsunami, J. Comput. Phys, 227, 7720–775,, 2008. 

Fernandez-Nieto, E., Parisot, M., Penel, Y., and Sainte-Marie, J.: A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows, Commun. Math. Sci., 16, 1169–1202, 2018. 

Fiedler, A. and Faleide, J. I.: Cenozoic sedimentation along the southwestern Barents Sea margin in relation to uplift and erosion of the shelf, Global Planet. Change, 12, 75–93,, 1996. 

Fine, I. V., Rabinovich, A. B., Bornhold, B. D., Thomson, R. E., and Kulikov, E. A.: The Grand Banks landslide-generated tsunami of November 18, 1929: preliminary analysis and numerical modelling, Mar. Geol., 215, 45–57,, 2005. 

Førland, E. J., Jacobsen, J. K. S., Denstadli, J. M., Lohmann, M., Hanssen-Bauer, I., Hygen, H. O., and Tømmervik, H.: Cool weather tourism under global warming: Comparing Arctic summer tourists' weather preferences with regional climate statistics and projections, Tourism Manage., 36, 567–579,, 2013. 

Fryer, G. J., Watts, P., and Pratson, L. F.: Source of the great tsunami of 1 April 1946: a landslide in the upper Aleutian forearc, Mar. Geol., 203, 201–218,, 2004. 

García, M., Ercilla, G., and Alonso, B.: Morphology and sedimentary systems in the Central Bransfield Basin, Antarctic Peninsula: sedimentary dynamics from shelf to basin, Basin Res., 21, 295–314,, 2009. 

García, M., Ercilla, G., Alonso, B., Casas, D., and Dowdeswell, J. A.: Sediment lithofacies, processes and sedimentary models in the Central Bransfield Basin, Antarctic Peninsula, since the Last Glacial Maximum, Mar. Geol., 290, 1–16,, 2011. 

González-Vida, J. M., Macías, J., Castro, M. J., Sánchez-Linares, C., de la Asunción, M., Ortega-Acosta, S., and Arcas, D.: The Lituya Bay landslide-generated mega-tsunami – numerical simulation and sensitivity analysis, Nat. Hazards Earth Syst. Sci., 19, 369–388,, 2019. 

Haflidason, H., Sejrup, H. P., Nygård, A., Mienert, J., Bryn, P., Lien, R., and Masson, D.: The Storegga Slide: architecture, geometry and slide development, Mar. Geol., 213, 201–234,, 2004. 

Haflidason, H., Lien, R., Sejrup, H. P., Forsberg, C. F., and Bryn, P.: The dating and morphometry of the Storegga Slide, Mar. Petrol. Geol., 22, 123–136,, 2005. 

Hampel, A., Hetzel, R., Maniatis, G., and Karow, T.: Three-dimensional numerical modelling of slip rate variations on normal and thrust fault arrays during ice cap growth and melting, J. Geophys. Res., 114, 406,, 2009. 

Harbitz, C. B., Løvholt, F., Pedersen, G., and Masson, D. G.: Mechanisms of tsunami generation by submarine landslides: a short review, N. J. Geol., 86, 255–264, 2006. 

Huggel, C., Zgraggen-Oswald, S., Haeberli, W., Kääb, A., Polkvoj, A., Galushkin, I., and Evans, S. G.: The 2002 rock/ice avalanche at Kolka/Karmadon, Russian Caucasus: assessment of extraordinary avalanche formation and mobility, and application of QuickBird satellite imagery, Nat. Hazards Earth Syst. Sci., 5, 173–187,, 2005. 

Iglesias, O., Lastras, G., Canals, M., Olabarrieta, M., Gonzalez, M., Aniel-Quiroga, Í., and De Mol, B.: The BIG'95 submarine landslide–generated tsunami: a numerical simulation, J. Geol., 120, 31–48,, 2012. 

Imamura, F., Boret, S. P., Suppasri, A., and Muhari, A.: Recent occurrences of serious tsunami damage and the future challenges of tsunami disaster risk reduction, Progress in Disaster Science, 1, 100009,, 2019. 

Jakobsson, M., Mayer, L., Coakley, B., Dowdeswell, J. A., Forbes, S., Fridman, B., Weatherall, P., Hodnesdal, H., Noormets, R., Pedersen, R., Rebesco, M., Scheke, H. W., Zarayskaya, Y., Accettella, D., Armstrong, Anderson, R. M., Bienhoff, P., Camerlenghi, A., Church, I., Edwards, M., Gardner, J. V., Hall, J. K., Hell, B., Hestvik, O., Krstoffersen, Y., Marcussen, C., Mohammad, R., Mosher, D., Nghiem, S. V., Pedrosa, M. T., Trabaglini, P. G., and Weatherall, P.: The International Bathymetric Chart of the Arctic Ocean (IBCAO) version 3.0, Geophys. Res. Lett., 39, L12609,, 2012. 

Knies, J., Matthiessen, J., Vogt, C., Laberg, J. S., Hjelstuen, B. O., Smelror, M., Larsen, E., Andreassen, K., Eidvin, T., Vorren, T. O.: The Plio-Pleistocene glaciation of the Barents Sea–Svalbard region: a new model based on revised chronostratigraphy, Quaterny Sci. Rev., 28, 812–829,, 2009. 

Kvalstad, T. J., Andresen, L., Forsberg, C. F., Berg, K., Bryn, P., and Wangen, M.: The Storegga slide: evaluation of triggering sources and slide mechanics, Marine and Petroleum Geology, Elsevier, 22, 245–256,, 2005. 

Laberg, J. S. and Vorren, T. O.: A Late Pleistocene submarine slide on the Bear Island Trough Mouth Fan, Geo-Mar. Lett., 13, 227–234,, 1993. 

Laberg, J. S., Vorren, T. O., and Knutsen, S. M.: The Lofoten contourite drift off Norway, Mar. Geol., 159, 1–6,, 1999. 

Laberg, J. S., Vorren, T.O., Dowdeswell, J. A., Kenyon, N. H., and Taylor, J.: The Andøya Slide and the Andøya Canyon, north-eastern Norwegian–Greenland Sea, Mar. Geol., 162, 259–275, 2000. 

Laberg, J. S., Vorren, T. O., Mienert, J., Evans, D., Lindberg, B., Ottesen, D., and Henriksen, S.: Late Quaternary palaeoenvironment and chronology in the Trænadjupet Slide area offshore Norway, Mar. Geol., 188, 35–60,, 2002. 

Lee, H. J.: Timing of occurrence of large submarine landslides on the Atlantic Ocean margin, Mar. Geol., 264, 53–64,, 2009. 

Lee, S., Wolberg, G., and Shin, S. Y.: Scattered data interpolation with multilevel B-splines, IEEE T. Vis. Comput. Gr., 3, 228–244,, 1997. 

L'Heureux, J. S., Vanneste, M., Rise, L., Brendryen, J., Forsberg, C. F., Nadim, F., Longva, O., Chand, S., Kvalstad, T. J., and Haflidason, H.: Stability, mobility and failure mechanism for landslides at the upper continental slope off Vesterålen, Norway, Mar. Geol., 346, 192–207,, 2013. 

León, R., Llorente, M., and Giménez-Moreno, C. J.: Marine Gas Hydrate Geohazard Assessment on the European Continental Margins. The Impact of Critical Knowledge Gaps, Appl. Sci., 11, 2865,, 2021. 

Lindberg, B., Laberg, J. S., and Vorren, T. O.: The Nyk Slide–morphology, progression, and age of a partly buried submarine slide offshore northern Norway, Mar. Geol., 213, 277–289,, 2004. 

Llopart, J.: Storfjorden Trough Mouth Fan (Western Barents Sea): Slope Failures in Polar Continental Margins; Significance of Stress Changes and Fluid Migration Induced by Glacial Cycles, PhD thesis, Universitat de Barcelona, 231 pp., (last access: 20 September 2022), 2016. 

Llopart, J., Urgeles, R., Camerlenghi, A., Lucchi, R. G., Rebesco, M., and De Mol, B.: Late Quaternary development of the Storfjorden and Kveithola Trough Mouth Fans, northwestern Barents Sea, Quaterny Sci. Rev. 129, 68–84,, 2015. 

Llopart, J., Urgeles, R., Forsberg, C. F., Camerlenghi, A., Vanneste, M., Rebesco, M., and Lantzsch, H.: Fluid flow and pore pressure development throughout the evolution of a trough mouth fan, western Barents Sea, Basin Res., 31, 487–513,, 2019. 

Løvholt, F., Pedersen, G., Harbitz, C. B., Glimsdal, S., and Kim, J.: On the characteristics of landslide tsunamis, Philos. T. R. Soc. A, 373, 2053,, 2015. 

Løvholt, F., Glimsdal, S., and Harbitz, C. B.: On the landslide tsunami uncertainty and hazard, Landslides, 17, 2301–2315,, 2020. 

Lucchi, R. G., Pedrosa, M. T., Camerlenghi, A., Urgeles, R., De Mol, B., and Rebesco, M.: Recent submarine landslides on the continental slope of Storfjorden and Kveithola Trough-Mouth Fans (North West Barents Sea), in: Submarine Mass Movement and Their Consequences, Advances in Natural and Technological Hazards Research, edited by: Yamada, Y., Kawamura, K., Ikehara, K., Ogawa, Y., Urgeles, R., Mosher, D., Chaytor, J., and Strasser, M., 31, Springer, Dordrecht, the Netherlands, 735–745,, 2012. 

Lucchi, R. G., Camerlenghi, A., Rebesco, M., Colmenero-Hidalgo, E., Sierro, F. J., Sagnotti, L., and Caburlotto, A.: Postglacial sedimentary processes on the Storfjorden and Kveithola trough mouth fans: Significance of extreme glacimarine sedimentation, Global Planet. Change, 111, 309–326,, 2013. 

Macías, J., Vázquez, J. T., Fernández-Salas, L. M., González-Vida, J. M., Bárcenas, P., Castro, M. J., Díaz-del-Río, V., and Alonso, B.: The Al-Borani submarine landslide and associated tsunami, A modelling approach, Mar. Geol, 361, 79–95,, 2015. 

Macías, J., Castro, M. J., Ortega, S., and González-Vida, J. M.: Performance assessment of Tsunami-HySEA model for NTHMP tsunami currents benchmarking, Field cases, Ocean Model., 152, 101645,, 2020. 

Macías, J., Escalante, C., and Castro, M. J.: Multilayer-HySEA model validation for landslide-generated tsunamis – Part 1: Rigid slides, Nat. Hazards Earth Syst. Sci., 21, 775–789,, 2021. 

Maslin, M., Mikkelsen, N., Vilela, C., and Haq, B.: Sea-level–and gas-hydrate–controlled catastrophic sediment failures of the Amazon Fan, Geology, 26, 1107–1110,<1107:SLAGHC>2.3.CO;2, 1998. 

Meleshko, V. P., Golitsyn, G. S., Govorkova, V. A., Demchenko, P. F., Eliseev, A. V., Kattsov, V. M., and Sporyshev, P. V.: Anthropogenic climate change in Russia in the 21st century: An ensemble of climate model projections, Meteorologiya i gidrologiya, 4, 38–49, 2004. 

Moernaut, J., Van Daele, M., Strasser, M., Clare, M. A., Heirman, K., Viel, M., Cardenas, J., Kilian, R., Ladrón de Guevara, B., Pino, M., Urrutia, R., and De Batist, M.: Lacustrine turbidites produced by surficial slope sediment remobilization: A mechanism for continuous and sensitive turbidite paleoseismic records, Mar. Geol, 384, 159–176,, 2017. 

Newton, A. M. and Huuse, M.: Glacial geomorphology of the central Barents Sea: implications for the dynamic deglaciation of the Barents Sea Ice Sheet, Mar. Geol, 387, 114–131,, 2017. 

Ottesen, D., Dowdeswell, J. A., and Rise, S.: Submarine landforms and the reconstruction of fast-flowing ice streams within a large Quaternary ice sheet: The 2500-km-long Norwegian-Svalbard margin (57–80 N), Geol. Soc. Am. Bull., 117, 1033–1050,, 2005. 

Pedrosa, M. T., Camerlenghi, A., De Mol, B., Urgeles, R., Rebesco, M., and Lucchi, R. G.: Seabed morphology and shallow sedimentary structure of the Storfjorden and Kveithola trough-mouth fans (North West Barents Sea), Mar. Geol., 286, 65–81,, 2011. 

Piper, D. J., Cochonat, P., and Morrison, M. L.: The sequence of events around the epicentre of the 1929 Grand Banks earthquake: initiation of debris flows and turbidity current inferred from sidescan sonar, Sedimentology, 46, 79–97,, 1999. 

Pirli, M., Schweitzer, J., and Paulsen, B.: The Storfjorden, Svalbard, 2008–2012 aftershock sequence: Seismotectonics in a polar environment, Tectonophysics, 601, 192–205,, 2013. 

Pouliquen, O. and Forterre, Y.: Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane, J. Fluid Mech., 453, 133–151,, 2002. 

Pudasaini, S. P.: Dynamics of submarine debris flow and tsunami, Acta Mech., 225, 2423–2434,, 2014. 

Ramadan, K. T.: Near- and far-field tsunami waves, displaced water volume, potential energy and velocity flow rates by a stochastic submarine earthquake source model, Nat. Hazards Earth Syst. Sci. Discuss. [preprint],, 2018. 

Rebesco, M., Pedrosa, M. T., Camerlenghi, A., Lucchi, R. G., Sauli, C., De Mol, B., Madrussani, G., Urgeles, R., Rossi, G., and Böhm, G.: One million years of climatic generated landslide events on the northwestern Barents Sea continental margin, in: Submarine Mass Movements and Their Consequences. Advances in Natural and Technological Hazards Research, edited by: Yamada, Y., Kawamura, K., Ikehara, K., Ogawa, Y., Urgeles, R., Mosher, D., Chaytor, J., and Strasser, M., 31, Springer, Dordrecht, the Netherlands, 747–756,, 2012. 

Rebesco, M., Laberg, J. S., Pedrosa, M. T., Camerlenghi, A., Lucchi, R. G., Zgur, F., and Wardell, N.: Onset and growth of trough-mouth fans on the north-western Barents Sea margin–implications for the evolution of the Barents Sea/Svalbard ice sheet, Quaternary Sci. Rev., 92, 227–234,, 2014. 

Rodríguez-Morata, C., Villacorta, S., Stoffel, M., and Ballesteros-Cánovas, J. A.: Assessing strategies to mitigate debris-flow risk in Abancay province, south-central Peruvian Andes, Geomorphology, 342, 127–139,, 2019. 

Salaree, A. and Okal, E.: Effects of bathymetry complexity on tsunami propagation: a spherical harmonics approach, Geophys. J. Int., 223, 632–647,, 2020. 

Sejrup, H. P., Hjelstuen, B. O., Torbjørn Dahlgren, K. I., Haflidason, H., Kuijpers, A., Nygård, A., Praeg, D., and Vorren, T. O.: Pleistocene glacial history of the NW European continental margin, Mar. Petrol. Geol., 22, 1111–1129,, 2005. 

Sierro, F. J., Andersen, N., Bassetti, M. A., Berné, S., Canals, M., Curtis, J. H., Dennielou, B., Flores, J. A., Frigola, J., Gonzalez-Mora, B., Grimalt, J. O., Hodell, D. A., Jouet, G., Pérez-Folgado, M., and Schneider, R.: Phase relationship between sea level and abrupt climate change, Quaterny Sci. Rev., 28, 2867–2881,, 2009. 

Skogseth, R., Olivier, L. L., Nilsen, F., Falck, E., Fraser, N., Tverberg, V., and Falk-Petersen, S.: Variability and decadal trends in the Isfjorden (Svalbard) ocean climate and circulation – An indicator for climate change in the European Arctic, Prog. Oceanogr., 187, 102394,, 2020. 

Solheim, A., Berg, K., Forsberg, C. F., and Bryn, P.: The Storegga Slide complex: repetitive large scale sliding with similar cause and development, Mar. Petrol. Geol., 22, 97–107,, 2005. 

Sun, Q. and Leslie, S.: Tsunamigenic potential of an incipient submarine slope failure in the northern South China Sea, Mar. Petrol. Geol., 112, 104111,, 2020. 

Synolakis, C. E., Yalciner, A. C., Borrero, J. C., and Plafker, G.: Modeling of the November 3, 1994 Skagway, Alaska tsunami, in: Solution to Coastal Disasters Conference, San Diego, California, United States, 24–27 February, 915–927,, 2002. 

Tappin, D. R.: Submarine mass failures as tsunami sources: their climate control, Philos. T. R. Soc. A, 368, 2417–2434,, 2010. 

Thomson, R. E., Rabinovich, A. B., Kulikov, E. A., Fine, I. V., and Bornhold, B. D.: On numerical simulation of the landslide-generated tsunami of November 3, 1994 in Skagway Harbor, Alaska, in: Tsunami research at the end of a critical decade, Springer, Dordrecht, 243–282,, 2001. 

Tinti, S., Bortolucci, E., and Romagnoli, C.: Computer simulations of tsunamis due to sector collapse at Stromboli, Italy, J. Volcanol. Geoth. Res., 96, 103–128,, 2000. 

Urgeles, R., Bahk, J. J., Lee, S. H., Horozal, S., Cukur, D., Kim, S. P., and Um, I. K.: Tsunami hazard from submarine landslides: scenario-based assessment in the Ulleung Basin, East Sea (Japan Sea), Geosci. J., 23, 439–460,, 2019. 

Urlaub, M., Talling, P. J., and Masson, D. G.: Timing and frequency of large submarine landslides: implications for understanding triggers and future geohazard, Quaterny Sci. Rev., 72, 63–82,, 2013. 

Vanneste, M., Mienert, J., and Bünz, S.: The Hinlopen Slide: a giant, submarine slope failure on the northern Svalbard margin, Arctic Ocean, Earth Planet. Sc. Lett., 245, 373–388,, 2006. 

Vanneste, M., Sultan, N., Garziglia, S., Forsberg, C. F., and L'Heureux, J. S.: Seafloor instabilities and sediment deformation processes: The need for integrated, multi-disciplinary investigations, Mar. Geol., 352, 183–214,, 2014. 

Vázquez, J. T., Ercilla, G., Alonso, B., Peláez, J. A., Palomino, D., León, R., Barcenas, P., Casas, D., Estrada, F., Fernández-Puga, M. C., Galindo-Zaldívar, J., Henares, J., Llorente, M., Sánchez-Guillamón, O., d'Acremont, E., Ammar, A., Chourak, M., Fernández-Salas, L. M., López-González, N., and Lafuerza, S.: Triggering Mechanisms of Tsunamis in the Gulf of Cadiz and the Alboran Sea: An Overview, in: Historical Earthquakes, Tsunamis and Archaeology in the Iberian Peninsula (Natural Science in Archaeology), edited by: Álvarez-Martí-Aguilar, M., and Prieto, F. M., Springer, Singapore, 65–104,, 2022.  

Winkelmann, D., Geissler, W., Schneider, J., and Stein, R.: Dynamics and timing of the Hinlopen/Yermak Megaslide north of Spitsbergen, Arctic Ocean, Mar. Geol,, 250, 34–50,, 2008. 

Worsley, D., Aga, O. J., Dalland, A., Elverhøi, A., and Thon, A.: The Geological History of Svalbard, Evolution of an Arctic Archipelago, Statoil, Asker Trykkeri, Stavanger, 121, 1986. 

Zhang, W., Almgren, A., Beckner, V., Bell, J., Blaschke, J., Chan, C., and Zingale, M.: AMReX: a framework for block-structured adaptive mesh refinement, J. Open Source Softw., 4, 1370–1370,, 2019. 

Short summary
The L-ML-HySEA (Landslide Multilayer Hyperbolic Systems and Efficient Algorithms) model of the tsunami triggered by the Storfjorden LS-1 landslide provides new insights into the sliding mechanism and bathymetry controlling the propagation, amplitude values and shoaling effects as well as coastal impact times. This case study provides new perspectives on tsunami hazard assessment in polar margins, where global climatic change and its related ocean warming may contribute to landslide trigger.
Final-revised paper