## Abstract

Sediment slumps are known to have generated important tsunamis such as the 1998 Papua New Guinea (PNG) and the 1929 Grand Banks events. Tsunami modellers commonly use solid blocks with short run-out distances to simulate these slumps. While such methods have the obvious advantage of being simple to use, they offer little or no insight into physical processes that drive the events. The importance of rotational slump motion to tsunamigenic potential is demonstrated in this study by employing a viscoplastic landslide model with Herschel–Bulkley rheology. A large number of simulations for different material properties and landslide configurations are carried out to link the slump's deformation, rheology, its translational and rotational kinematics, to its tsunami genesis. The yield strength of the slump is shown to be the primary material property that determines the tsunami genesis. This viscoplastic model is further employed to simulate the 1929 Grand Banks tsunami using updated geological source information. The results of this case study suggest that the viscoplastic model can be used to simulate complex slump-induced tsunami. The simulations of the 1929 Grand Banks event also indicate that a pure slump mechanism is more tsunamigenic than a corresponding translational landslide mechanism.

Landslides constitute the second-most important tsunami source worldwide after earthquakes (Tappin 2010; Harbitz *et al.* 2014; Yavari-Ramshe and Ataie-Ashtiani 2016). Most recently, the 2018 Anak Krakatoa event caused several hundred fatalities (Grilli *et al.* 2019). Between 2007 and 2017 a string of at least five additional large subaerial landslides impacted water and generated run-up heights in the range of 30 to 150 m (Sepúlveda and Serey 2009; Wang *et al.* 2015; George *et al.* 2017; Gylfadóttir *et al.* 2017; Paris *et al.* 2019). Submarine landslide tsunamis are less frequent than these subaerial landslide tsunamis, but the largest recognized events worldwide indisputably illustrate their destructive potential and importance for society. Fatal examples of such submarine landslides are the 1998 Papua New Guinea (PNG) (Synolakis *et al.* 2002), 1992 Flores Island (Yeh *et al.* 1993), 1979 Lembata Island (Yudhicara *et al.* 2015) and 1929 Grand Banks landslides (Løvholt *et al.* 2019).

Slumps constitute a subset of landslides that are typically characterized by a rotational impulsive slope failure, a relatively coherent mass displacement, and a short landslide run-out distance. At least two of the above-mentioned events, the 1998 PNG and the 1929 Grand Banks events, were caused by rotational slumps. The study of the PNG event also led to acknowledgement in the scientific community that submarine slumps can cause large tsunamis (Bardet *et al.* 2003; Tappin *et al.* 2008). This tsunami has been successfully modelled using an approach where the landslide motion is a rigid block that follows a prescribed motion (Synolakis *et al.* 2002; Tappin *et al.* 2008), by tuning the block motion to comply with wave observations. A similar approach was adopted for modelling the slump part of the 1929 Grand Banks event (Løvholt *et al.* 2019). The rigid block approach was successful in these studies, because the block could mimic the rotational motion of the slump causing the tsunami genesis in an idealized and simple way, but did not include the updated geological source information from Schulten *et al.* (2019*b*), which envisaged a slump that partly evacuated the source area. Although this block modelling approach can help to shed light on the slump motion of past events, it has several obvious shortcomings. Firstly, this method does not include landslide deformation effects that are evident from geophysical data. Secondly, these models cannot be used to take into account the landslide material properties such as the yield strength, and its effect on the landslide dynamics.

Recent modelling efforts show that the landslide rheology and deformation is important for quantifying and understanding landslide tsunami genesis (Løvholt *et al.* 2017; Yavari-Ramshe and Ataie-Ashtiani 2019). Traditionally, such landslide tsunami studies are based on translational landslide models. However, translational landslides are believed to give rise to a different generation mechanism than slumps, as they do not exhibit a rotational motion as slumps do (Løvholt *et al.* 2015). Until recently, slump models that include a more sophisticated deformation and rheology had not been applied to slump-induced tsunamis. Schambach *et al.* (2018) provided back-to-back analysis with a viscous landslide model and a rigid block model simulating slumps, with both models showing similar results. Ren *et al.* (2019) used a viscoplastic landslide model to generate the slump tsunami due to the 1998 PNG failure, with simulation results that compare favourably with tsunami inundation observations. These studies (Schambach *et al.* 2018; Ren *et al.* 2019) show that a slump tsunami can be effectively modelled using a landslide dynamics model. This method allows for a more flexible, general modelling treatment of the slump tsunami genesis, including material properties, deformation and complex topography, which will be utilized herein.

In this paper, we will use the viscoplastic model BingClaw (Løvholt *et al.* 2017; Kim *et al.* 2019), coupled to the dispersive long-wave solver GloBouss (Løvholt *et al.* 2008), to study slump-induced tsunamis. We will first study landslide dynamics and tsunami genesis in an idealized geometry in one-horizontal dimension (1HD). The main aims of this idealized study are, for the first time, to:

quantify relationships between landslide material yield strength, the resulting slump kinematics and dynamics, and slump tsunamigenic potential; and

identify the extent to which slump tsunamigenic potential can be attributed to translational and rotational slump kinematics, such as the angular momentum.

We will apply the same model setup in two horizontal dimensions (2HD) to study a real case, namely the 1929 Grand Banks landslide and tsunami. The main emphasis of the real case example is to ensure that the landslide parameters and settings in the idealized study can yield a realistic range of analysis. However, a detailed study of the event is left for future investigations.

## The 1929 Grand Banks landslide and tsunami

On 8 November 1929 a M_{w} 7.2 earthquake caused a massive landslide on the Grand Banks south of Newfoundland (Heezen and Ewing 1952; Piper *et al.* 1999) (see Fig. 1). This submarine mass failure comprises by far the largest landslide volume (*c.* 500 km^{3}) in historical time, worldwide. Deposits far from the landslide failure area and cable breaks (Heezen *et al.* 1954) suggest that the landslide evolved into a turbidity current. The landslide caused a tsunami several metres high at the Burin Peninsula on the south coast of Newfoundland, and waves were also recorded along the entire US East Coast, Bermuda and the Azores (Fine *et al.* 2005). Initial field evidence of the landslide deposits suggested that only turbidity current masses were available in the far field (Schulten *et al.* 2019*a*). Piper *et al.* (1999) noted that the Grand Banks landslide was a widely distributed surficial sediment failure, and Mosher and Piper (2007) noted from newly acquired multibeam bathymetric data that there was no evidence of a massive slump failure on the St Pierre Slope. As the turbidity current itself is likely not the cause of the tsunami, it has been difficult to link the tsunami genesis directly to landslide field evidence. Based on new field investigations of the slope failure, however, Schulten *et al.* (2019*a*) and Løvholt *et al.* (2019) suggested that the near-field tsunami was caused by a massive slump. Løvholt *et al.* (2019) further hypothesized that the more widespread near-surface landslide failure as mapped by Piper *et al.* (1999) and Schulten *et al.* (2019*a*) caused the far-field tsunamis, and that the landslide possibly disintegrated into the turbidity current. Løvholt *et al.* (2019) used a simplified block source and a slump volume of 17 km^{3} to model the slump. However, the analysis of newly identified faults and horizons in the St Pierre Slope by Schulten *et al.* (2019*b*) suggest a much larger slump volume of *c.* 390 km^{3} for the primary southward slump motion. This new interpretation for the 1929 Grand Banks slump is crucial for testing whether or not our viscoplastic flow model is suited to simulate slumps. Moreover, Schulten *et al.* (2019*b*) suggest that the slump was not confined only between the structural faults containing the slump mass, but also that parts of the landslide transgressed the downslope end of the slump source area through the channel systems, which is different from the assumption of Schulten *et al.* (2019*a*) and Løvholt *et al.* (2019).

## Methods

### Landslide model

In this paper, the viscoplastic landslide model BingClaw (Løvholt *et al.* 2017; Kim *et al.* 2019; Vanneste *et al.* 2019) is used to simulate the slump dynamics. The model implements the Herschel–Bulkley rheology in a two-layer depth-averaged formulation. Under simple shear conditions, the shear strain in the Herschel–Bulkley fluid is described as:
*μ*. *τ* and *τ*_{y} are shear stress and yield strength, respectively, and *n* the flow exponent. For a detailed description and derivation of the model, see Kim *et al.* (2019).

BingClaw solves the mass conservation equation integrated over the landslide depth (equation 3), the momentum conservation equation integrated separately over the plug layer depth (equation 4), and shear layer depth (equation 5), in 2HD. The unknown variables are bed-normal plug layer thickness *d*_{p}, bed-normal shear layer thickness *d*_{s}, plug layer volume flux per unit length *d* = *d*_{p} + *d*_{s} is the total thickness of the layers. Indices ‘p’ and ‘s’ indicate plug and shear layer, respectively (see Fig. 2).
*C*_{m} is the added-mass coefficient, *ρ*_{w} the density of ambient water, *ρ*_{d} the density of the slump material, *α* the velocity form factor, and *t* the time coordinate. The reduced gravitational acceleration is given by *g* is the gravitational acceleration, *b* is the bathymetric depth, *τ*_{d} is the viscous drag at the free surface, split into a skin friction term *τ*_{f} given by
*τ*_{p} given by
*C*_{F} and *C*_{P} are skin friction and pressure drag coefficients, respectively, and the viscous contribution of the net shear stress at the bed is given by *τ*_{y}*f*_{s} where
*β* is a shape factor depending on the rheological flow exponent *n* (Huang and Garcia 1998; Imran *et al.* 2001; Kim *et al.* 2019).

BingClaw combines a finite volume method for the leading-order terms with a finite difference model for the source terms. The model is implemented employing the conservation law package ClawPack (Mandli *et al.* 2016) using the GeoClaw module (Berger *et al.* 2011). If the earth pressure *p* = *ρ*_{d}*g*′*d*∇(*d* + *b*) does not exceed the material's shear strength in a given computational cell, no motion is imposed in that cell. Otherwise a Godunov fractional step method is used for the dynamic equations. First the equations without friction terms are solved using the finite volume method in ClawPack, then the frictional terms are accounted for the next fractional step.

### Tsunami model

We use the dispersive long wave model GloBouss (Løvholt *et al.* 2008, 2010; Pedersen and Løvholt 2008) to propagate the tsunami over varying bathymetry. In this study, we only use the model in linearized mode as we mainly study the tsunami in deep water, where non-linearities are unimportant.

When terms and factors that are not used herein are omitted (non-linear terms, Coriolis terms, spherical coordinate map-factors and dispersion enhancement terms), the hydrodynamic equations used in this paper read
*q* is a source flux term, which relates the landslide model to the tsunami model through the landslide volumetric displacement (explained below); *h* is the water depth relative to the mean sea-surface elevation, *η* the sea-surface elevation and

In GloBouss the equations are discretized on a staggered C-grid (Mesinger and Arakawa 1976) in space and time to give an implicit finite difference method. An alternating direction implicit method (ADI) is used for solving the implicit algebraic equation systems for each time step. The model does not incorporate features like drying or wetting, so we cannot use this model to simulate dry-land inundation.

The slump causes a temporal volumetric change of the bathymetry, which is the primary source for the tsunami genesis. These source fields are then run through a low pass filter that conveys seabed displacements to sea-surface displacements based on full potential wave theory (Kajiura 1963; Løvholt *et al.* 2015) that transfers ∂*d*/∂*t* into *q*(*x*, *y*, *t*).

### Model setup

The geometrical setup is based on the most recent 1929 Grand Banks landslide information provided by Schulten *et al.* (2019*b*). Our first objective is to link a rotational slump motion to tsunami genesis in a systematic fashion, where the slump is confined between an upslope and a downslope fault. To force the slump to stay between these structures, we choose to excavate the slump mass from the seabed, replace it with our viscoplastic material for the initial setup, and elevate the face of the downslope fault (see Fig. 3a, b). While we acknowledge that this geometrical description would likely differ from more complex field observations, this was a necessary simplification to force the viscoplastic material not to evacuate the structure. To this end, we first simulate the slump tsunami in 1HD. The aim is to study idealized effects of kinematics and landslide parameters on tsunami genesis. Secondly, we study a 2HD scenario for the 1929 Grand Banks event, for which the purpose is to provide a realistic parameter range for the 1HD study.

#### 1HD study

The 1HD geometries applied here are simplified from slope transects taken from the general Laurentian Fan bathymetry. As shown in Figure 3a and b, different bathymetries are investigated to study the sensitivity to the slope configuration of the slump source. The bathymetry outside the slump towards the shore is gentler with a constant inclination of 0.05° in all cases.

The computational domain for the landslide model has a total length of 50 km in the *x*-direction with a spatial resolution of Δ*x* = 80 m. However, due to the computational stencil of BingClaw, several cells in the azimuthal *y*-direction are required. Non-reflecting outflow conditions are applied at the boundaries. The Kajiura-type full potential filter is run over the same length as the landslide model. Grid resolution for the Kajiura filtered output is also 80 m. For GloBouss we cover a computational domain extending 450 km horizontally, and with a resolution of 220 m. We apply a sponge layer at the right boundary, from 250 to 450 km, that relaxes the offshore-going waves (Pedersen and Løvholt 2008). No-flux conditions are applied at the other boundaries. In GloBouss a 1HD computation involves a single wet row of cells between two dry rows of ghost cells. Spatial and temporal grid refinement tests on the landslide model BingClaw, the full potential Kajiura-type filter, and the tsunami model GloBouss are described in Appendix A.

Default model input parameters (i.e. density and hydrodynamic resistance) are listed in Table 1, and geometrical and geotechnical model input parameters used for the sensitivity analysis are listed in Table 2. In order to establish the list of landslide input parameters, we ran several simulations to achieve a full parameter range that spans the relevant sensitivity range for tsunami genesis. By combining all relevant geotechnical and hydrodynamic resistance parameters for each geometrical setting, we ended up running 2640 simulations for the sensitivity analysis. We refer to 1440 model runs with constant slump volumes per unit width and variable initial slump surface slope angle as set S1, and 1440 model runs with a constant initial slump surface slope angle and variable volumes per unit width as set S2; 240 simulations overlap in set S1 and set S2.

A simplified basic geometry is defined by an initial slope angle *θ* = 2.5°, as retrieved from the Laurentian fan, and a volume per unit width *A* = 5.2 km^{2}, which multiplied with a slump width of *W* = 33 km yields a total volume *V* = 175 km^{3} as suggested by Schulten *et al.* (2019*b*) for the upper part of the 1929 Grand Banks slump. Then, in simulation set S1, *θ* is varied between 1 and 3.5°, while keeping *A* constant. Likewise, in set S2, *A* is varied between 1.7 km^{2} and 7.5 km^{2}, while keeping *θ* constant. In each case the parabolic shape of the rigid seabed is adjusted accordingly.

For a very soft slump material (e.g. low values of *τ*_{y} in Table 2), the mass can be so mobile that it artificially reflects from the lower fault face and propagates back upslope and may even continue to slosh back and forth. This spurious sloshing occurs partly due to simplifications in the applied slump model, partly due to the geometrical setup and partly due to too small values employed for the landslide strength. Time-series of two examples of the centre-of-mass motions, which is used to filter events, are shown in Figure 4. The centre-of-mass velocities have a smoother time evolution than the maximum velocities. If an event gets a negative centre-of-mass velocity, it is removed from the analysis to avoid the artificial sloshing. This criterion was based on analyses of the wave generation for the sloshing events, where it transpired that events with negative centre-of-mass velocities influenced the wave generation significantly. An example of the artificial scaling behaviour that can be expected is discussed in one of our analyses below. The number of non-sloshing events as well as events where the yield strength is too large for the mass to mobilize the landslide (i.e. stable sediments), are shown for both set S1 and set S2 in Figure 5.

#### 2HD study

The slump configuration with the new information provided by Schulten *et al.* (2019*b*) is used to simulate the slump dynamics. We distinguish between two different scenarios, an over-topping (where a part of the material escapes in the lower extremity) and a pure slump. For the pure slump the mass is confined to a source area limited by walls at the downslope extremity and at the two sides. It generates a rotational slump motion in a similar way as in the 1HD study (see Fig. 3c). We note that an over-topping scenario is considered as most likely by Schulten *et al.* (2019*b*) (see Fig. 3d). In the case of over-topping, the model geometry is set up to allow the slide material to continue as a translational landslide outside the region of mass failure. The further disintegration into the turbidity current observed in the field is, however, not included in the model. We note that the main orientation of this slump geometry is southward, which was also assumed by Løvholt *et al.* (2019). Yet, the revised slump volume used in the 2HD analysis here (390 km^{3}) is considerably larger than what was assumed by Løvholt *et al.* (2019).

Model bathymetries are based on the online geographical GEBCO 2014 Grid with 463 m cell size in longitude and latitude. The depth matrix for the landslide and source computations covers a rectangle with lengths of 114 and 255 km in the longitudinal and latitudinal directions, respectively. For the landslide model, a grid resolution of 185 m is used, while a resolution of 463 m suffices for the surface response. As in the 1HD slump model, there is non-reflecting outflow at the four boundaries. The grid for the tsunami computations is larger and covers a rectangle of 616 km (longitude) by 555 km (latitude). It has a resolution of 463 m and includes the source area, the southern coast of Newfoundland and the eastern coast of Nova Scotia (see Fig. 1). At all four boundaries we apply a sponge layer of 22 km width where the waves are relaxed and apply a minimum computation depth of 10 m in order to avoid spurious oscillations in shallow waters. Spatial and temporal grid refinement tests on the landslide model BingClaw, the full potential Kajiura-type filter and the tsunami model GloBouss are discussed in Appendix A. Default model parameters are presented in Table 1, geotechnical and geometrical parameters are given in bold in Table 2.

## Results

### 1HD parametric sensitivity study

#### Example of tsunami-genesis mechanism

We first analyse, through one single simulation, the slump tsunami-genesis mechanism. We use the following BingClaw parameters, namely *τ*_{y} = 70 kPa, *μ* = 10 kPa s* ^{n}*, and

*n*= 0.25. The slump surface slope angle is

*θ*= 2.5°, the slump volume per unit width is

*A*= 5.2 km

^{2}and the water depth of the initial centre of mass is

*c.*1750 m. At

*c.*236 km from shore, a maximum vertical landslide displacement of

*c.*100 m is obtained, which is similar to what was suggested by Schulten

*et al.*(2019

*b*). Figure 6a shows the slump motion at different times. The corresponding generated waves at different times are displayed in Figure 6b. While the slump mass rotates around its mass centre, the downslope part of the rotational slump pushes water upwards creating a positive wave at the surface, whereas the upslope part of the slump pulls down water and causes a trough at the surface.

Next, we ran two separate simulations for the same example, one using only the positive flux part of the slump source term, and another only using the negative flux part. Figure 7 shows the generated total wave (in solid lines), as well as the wave component only due to the slump uplift (in long dashed lines), and the wave component due to the slump depression (in short dashed lines). Both the generated wave elevations and wave troughs continuously split into landward and offshore travelling waves as long as the slump motion continues and add to the already propagated waves. Because the slump's upward and downward motions are spatially shifted, the landward wave-elevation travels slightly behind the landward wave-trough. Only a partial overlap of this wave-trough and wave-elevation occurs, which results in a landward trough followed by an elevation. The positive and negative amplitudes of this total wave, when travelled out of the source area, are roughly half of the maximum/minimum elevations from pure positive and negative source components.

This mechanism was discussed by Løvholt *et al.* (2005, 2015) and Haugen *et al.* (2005), but mainly for translational landslides. Based on analyses of the 1998 PNG event, Løvholt *et al.* (2015) suggested that the interaction between rear and frontal waves was limited for slumps, and that their wave generation was more efficient than for translational landslides. However, the present analysis shows that the interaction between the frontal and rear wave clearly reduces the maximum elevation of the total wave for the 1929 Grand Banks slump. We stress that for other slump configurations and material parameters the picture could be different.

#### Relationship between geotechnical parameters and tsunami genesis

Figure 8 shows the sensitivity of the maximum landward sea-surface elevation *η*_{max} to various input parameters. *η*_{max} is evaluated 900 s after the slump mass release such that the wave with the highest crest has propagated out of the source area. The various input parameters include the slump material's yield strength *τ*_{y}, the volume per unit width *A*, the initial slump surface slope angle *θ*, the dynamic consistency *μ*, and the flow exponent *n*. In all cases, *η*_{max} is plotted as a function of *τ*_{y}, and increases consistently with decreasing *τ*_{y}. As expected, *η*_{max} also increases with *θ* and *A*. Furthermore, we see that *η*_{max} is only moderately dependent on *μ*. The flow exponent *n* has a negligible influence on tsunami genesis, except when very small.

#### Relationship between landslide translational kinematics and tsunami genesis

Figure 9 shows relationships between maximum bed-parallel and vertical slump kinematics, and maximum and absolute minimum landward sea-surface elevations *η*_{max} and *η*_{min} for set S1. We recall that for S1, the initial slump surface slope angle *θ* is variable and the volume per unit width is constant at *A* = 5.2 km^{2}. The maximum kinematic quantities are calculated over the full computational domain for all times, whereas *η*_{max} and *η*_{min} are evaluated at a time of 900 s. Figure 9a shows scaled *η*_{max} and *η*_{min} as a function of the scaled maximum bed-parallel velocity *η*_{max} increases with *et al.* (2015) for slumps with small Froude numbers. The linear scaling relation should exist when there is no interaction between the frontal-wave elevation and rear-wave trough. However, in this case, there is clearly a destructive interference (see Fig. 7), which leads to a less effective wave generation. Figure 9c shows the relationship between the scaled maximum vertical velocity *η*_{max} and *η*_{min}. Unlike in Figure 9a, we do not observe a simple power-law relationship. There is also clearly more scatter in the vertical velocity plot. Further, processing of the kinematic output also verifies that maximum velocities, *η*_{max} shows a similar power-law dependency on *et al.* 2005, 2015). These studies concluded that the horizontal acceleration strongly influences tsunami genesis and, in particular, Løvholt *et al.* (2005, 2015) suggest the same linear relationship between *η*_{max} and

We recall that set S2 has a constant initial slump surface slope angle *θ* = 2.5°, but has different values for the volume per unit width *A*. The velocity is multiplied by the slump's total mass per width to quantify the momentum and to analyse how the momentum correlates with *η*_{max} and *η*_{min}. Figure 10a shows that *η*_{max} and *η*_{min} as functions of *η*_{max} is 0.9. Figure 10c shows that *η*_{max} and *η*_{min} have a similar relationship with the vertical maximum momentum *η*_{max} and *η*_{min} follow similar relationships as the ones derived for *η*_{max} and the rate of *η*_{max} and *η*_{min} against vertical momentum and momentum rates show less variability than the corresponding plots for *η*_{max} against vertical velocities and accelerations for set S1.

A Froude number, *Fr*, is defined as the maximum horizontal central mass velocity divided by the linear wave speed *H* = 2000 m. A nearly unitary *Fr* means the slump's horizontal central mass speed and the tsunami speed are the same, which represents the most efficient tsunami-genesis mechanism (Løvholt *et al.* 2015). In our study, *Fr* is invariably much smaller than unity. Figure 11 shows scaled *η*_{max} and *η*_{min} as a function of *Fr* for set S1, which represents the left side of the height–velocity curve peak in Ward (2001, fig. 3). We see that the growth rate of *η*_{max} as a function of *Fr* is slower than when we use the maximum landslide velocity (i.e. in Fig. 9a). On the other hand, we visually observe a slight misfit for the largest values of *Fr*, which may suggest that the exponent is not linear, possibly increasing with larger *Fr*. We note that Figure 11b also shows the results for the unfiltered simulations (i.e. including spurious sloshing events). Investigating Figure 11a and b, we see that the filter removes scenarios above *Fr* ≈ 0.13. For larger Froude numbers, the scaling of the unfiltered maximum landward sea-surface elevation *η*_{max} separates from the scaling of the absolute minimum sea-surface elevation *η*_{min}, and the separation occurs above *Fr* ≈ 0.15, say. The more rapid increase in the *η*_{max} with *Fr* is interpreted as a spurious result of the model (and hence filtered). On the other hand, we see that the scaling relationship for *η*_{min} is virtually unchanged for high Froude numbers (filtered events). The leading landward troughs are unaffected by the sloshing, which hints that a linear Froude scaling should also be expected for somewhat larger Froude numbers than those analysed elsewhere in this paper.

#### Relationship between landslide rotational kinematics and tsunami genesis

Slumps are mainly rotational and display different kinematics compared to translational landslides with long run-out. Here, we analyse to which extent the slump's scaled maximum angular momentum *L*_{max} is attributed to the slump's tsunamigenic potential. The technical derivation of this quantity is given in Appendix A. Figure 12a shows a power-law relationship between *L*_{max}, *η*_{max} and *η*_{min} for set S1. The exponent for *η*_{max} is 0.76. Figure 12b shows that the dependency between *L*_{max}, *η*_{max} and *η*_{min} for set S2 has significantly more scatter, and a less clear correlation. The fitted exponent is 0.66 for *η*_{max}. In both cases, the data exhibit little scatter for large *L*_{max}.

### 2HD study related to the 1929 Grand Banks event

#### Slump scenarios with over-topping

Figure 13 shows the simulated motion of the slump with over-topping for a volume of *V* = 390 km^{3} and a yield strength of *τ*_{y} = 85 kPa. At 300 s, the slump is still confined in the fault structure. Around *t* = 600 s the slump has its maximum vertical uplift of *c.* 400 m at its downslope extremity while parts of the slump mass escape the faulted pit and continue downslope as a translational landslide. This over-topping results in a 100 m high frontal landslide height. The output at 1380 s shows the landslide flowing into the Laurentian Fan region.

In the early phase, the generated wave (see Fig. 14) has a positive sea-surface elevation at the southern end of the slump area and a negative elevation at the northern end of the slump area. It is aligned north–south along the failure surface slope orientation. One hour after the slump mass release, the wave has started to turn gradually northwards and reaches the latitude 46° N after two hours. The main wave direction is towards the Burin Peninsula, whereas there is also a focus towards the Avalon Peninsula further east. Results extracted over the transect just south of Burin further show that maximum offshore sea-surface elevations range from 4 to 9 m for different landslide yield strengths (see Fig. 15a), which are in the same range or, for the lowest yield strengths, somewhat higher than those found by Løvholt *et al.* (2019). Figure 16a shows the maximum sea-surface elevations over the full simulation time, which coincides with the large waves observed near the Burin Peninsula (see, for example, Fine *et al.* 2005). Field observations of run-up elsewhere were mostly below 2 m; however, our simulations show waves just as large near Nova Scotia and the Avalon Peninsula as found near Burin.

More tuning would be necessary to provide a closer agreement with the data. For instance, Schulten *et al.* (2019*b*) found a vertical uplift of the slump mass of 100 m at its downslope extremity, although our example with *τ*_{y} = 85 kPa produces a much larger vertical uplift of *c.* 400 m. Our simulations merely provide a first attempt. However, the simulations clearly show that the viscoplastic model is capable of producing sufficiently strong slump-induced waves to produce a tsunami at least of the size of the 1929 Grand Banks event. We re-emphasize that our objective here was primarily to investigate whether the material parameter ranges for the 1HD case were representative of a real example, and this analysis shows that they are.

#### Slump scenarios without over-topping

We turn our attention to the pure slump, which is confined by the outreaching fault at the lower extremity, and to its tsunami using the same volume and material parameters as for the over-topping slump. Figure 17 shows the slump thickness 0, 300, 600 and 840 s after the mass release. At the last time the slump motion has stopped. The maximum vertical uplift is *c.* 800 m, which is twice as much as for the over-topping slump due to confinement. The spreading waves (see Fig. 18) and the total wave field (see Fig. 16b) have a similar radiation pattern as the over-topping slump tsunami; however, the positive generated waves are significantly larger along the ridges between the Laurentian Channel, the Halibut Channel and the Haddock Channel than the waves for the over-topping slump source (see Fig. 16). Figure 15b shows a transect just south of Burin with maximum offshore sea-surface elevations ranging from 1 to 5 m, which are, however, in the same range as the sea-surface elevations for the over-topping case. Near longitude 55.7°, we see that the over-topping scenario produces slightly larger waves than the pure slump. Still, on an overall basis, we suggest that the pure slump event seems to be a slightly more efficient tsunami generator than the over-topping event. This was confirmed by our own preliminary work on simulating Grand Banks (results not shown) with other slump configurations, where the difference was even clearer.

## Concluding remarks

In this paper, we have conducted a study of slump-induced tsunamis using a depth-averaged viscoplastic landslide model as the tsunami source, and a linear dispersive long-wave model for the tsunami propagation. Our main emphasis has been to study the sensitivity to slump material properties in 1HD on idealized geometries and the resulting slump kinematics on tsunami genesis. Contrary to most previous studies, our use of a viscoplastic landslide model allows us to link the tsunami directly to slump material properties, and avoid *ad-hoc* assumptions commonly made using a block model approach where the slump motion is prescribed. This refined model allows a more generalized treatment of slump sources, and hence is not limited to models that retrofit block source properties to simulate past events.

This study has shown that the material parameter that influences tsunami genesis the most is the initial yield strength of the sediment. Similar conclusions were reached for translational landslides in studies of the tsunami genesis of the Storegga landslide, for example (Kim *et al.* 2019). Moreover, our range of the dynamic landslide consistency (related to the viscosity) shows a more moderate influence on tsunami genesis. Naturally, geometrical factors, such as the slope angle and volume of the slump, were found to have a strong influence on the tsunami genesis too. Several kinematic properties were found to correlate well with the maximum landward sea-surface elevation. For the case of constant slide volume, the maximum landward sea-surface elevation increases monotonically with both scaled bed-parallel maximum velocity and acceleration mimicking a power-law relationship. The maximum landward sea-surface elevation also increases monotonically with vertical acceleration and velocity, but a less systematic relationship was found in this case. For the more general cases where variable volumes were investigated, the maximum bed-parallel momentum and momentum rates correlate well with the maximum landward sea-surface elevation, while the maximum landward sea-surface elevation had a somewhat less systematic relationship with corresponding vertical momentum and momentum rates.

Some of the findings of this study have been identified already in past studies (Tinti *et al.* 2001; Ward 2001; Løvholt *et al.* 2005), but only for translational landslides with a simplified block source representation. Here, we show that similar relationships between landslide velocities, accelerations and momentum apply for slumps. In particular, we find the scaling between the maximum height of the generated wave and the maximum bed-parallel landslide speed divided by the wave celerity, *η*_{max}∝*Fr*^{0.9}. We note that the exponent of 0.9 is less than the linear relationship (i.e. exponent 1) expected for small Froude numbers for frontal wave elevations and rear wave troughs without any interference (Løvholt *et al.* 2015). In our study, we clearly have destructive interference between the waves caused at the front and rear part of the slump which reduced the tsunamigenic potential. However, we find, similar to Løvholt *et al.* (2005, 2015), an almost linear scaling with the horizontal landslide acceleration, which is hence clearly a good proxy for tsunamigenic potential. An additional finding from our study is that the angular momentum shows a particularly good correlation with the maximum landward sea-surface elevation. This suggests that the tsunamigenic potential can be directly linked to rotational kinematic properties of the slump. We are unaware of previous studies that identify such a relationship.

A second part of the study is devoted to studying the 1929 Grand Banks slump and tsunami in a real topographical setting. This was primarily done to investigate whether the parameter ranges used in the viscoplastic slump model in 1HD were realistic. A detailed analysis of the 1929 Grand Banks event with an emphasis on obtaining a close match with field observations of the tsunami was not attempted. Nevertheless, our model was set up with new field observations by Schulten *et al.* (2019*b*) to illustrate how the geological interpretation provided a significantly revised explanation for the slump event. Schulten *et al.* (2019*b*) concluded that the 1929 Grand Banks slump failed mainly southwards, and that the main slump volume was much larger than previously thought (390 km^{3}). Our tsunami modelling suggests that a viscoplastic model indeed should be capable of producing sufficiently large waves. The 1929 Grand Bank event also served the purpose of testing how a complex event with slump failure and over-topping compares with a pure slump event with respect to tsunami genesis. We found that the pure slump produced larger overall waves compared to the over-topping scenario. All in all, the 1929 Grand Banks model including new field observations for the slump event and an idealized study in 1HD could revise our understanding of tsunami genesis.

## Acknowledgements

Computations were done on a computer located at the Norwegian Geotechnical Institute in Oslo, Norway. The authors are very grateful for positive criticism from the editor Aggeliki Georgiopoulou, and the reviewers David Mosher and Alberto Armigliato.

## Funding

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No.721403.

## Author contributions

**TZ**: data curation (lead), formal analysis (lead), investigation (lead), methodology (equal), software (equal), validation (equal), visualization (lead), writing – original draft (lead), writing – review & editing (equal); **FL**: conceptualization (lead), data curation (supporting), formal analysis (supporting), investigation (supporting), methodology (equal), project administration (lead), software (supporting), supervision (lead), validation (equal), writing – original draft (supporting), writing – review & editing (lead); **GKP**: formal analysis (supporting), investigation (supporting), methodology (equal), software (equal), supervision (equal), validation (equal), writing – original draft (supporting); **CBH**: conceptualization (supporting), data curation (supporting), formal analysis (supporting), investigation (supporting), methodology (equal), project administration (equal), supervision (equal), validation (equal), writing – review & editing (equal).

## Appendix A

### Grid refinement tests

For the 1HD simulations, we conducted grid refinement tests on the spatial grid for the slump model, tsunami model, and the Kajiura-type filter (resolutions and parameters in Table A1). For the slump model, we tested soft slump materials, low *τ*_{y} and low *µ*. The slump thickness depended strongly on the grid resolution for Δ*x* > 80 m. After 240 s, for instance, the slump thickness at the lower extremity is 8% thinner for Δ*x =* 160 m than for Δ*x* = 26 m, and for resolutions Δ*x* ⩽ 80 m, the slump thickness varies maximum by 4%. The slump thickness at the upslope part coincides for resolutions Δ*x* ⩽ 80 m, but gives twice the corresponding slump thickness for Δ*x =* 160 m. Thus, a spatial grid resolution of 80 m is chosen for further use. The time step, Δ*t*, is adapted during the simulation to keep the Courant–Friedrichs–Lewy (CFL) number (Courant *et al.* 1967),
*U*_{0} is the maximum particle speed in the slide body. In all our 1HD model runs, we use a CFL = 0.45, which yields stable behaviour (greatest landslide velocities are *c.* 70 m s^{−1}). When the source input is fed into the tsunami model each 30 s, we have a deviation of less than 2% from the smallest interval tested (5 s) at *t* = 480 s. Hence, we stay with 30 s. Since the surface response is smoother than the slide surface, application of the same spatial grid resolution for the Kajiura-type filter as for the landslide model, 80 m, is more than adequate. We tested spatial grid resolutions for the tsunami model 980 s after slump mass release. The maximum landward sea-surface elevation of a resolution of 220 m only deviated by 0.6% from the elevation of a finer resolution of 55 m. Therefore, we further used the 220 m resolution. The CFL number used is 0.5.

In 2HD, we executed spatial and temporal grid refinement tests of the landslide model BingClaw, the Kajiura-type filter, and the tsunami model GloBouss. All numerical parameters can be found in Table A2. For BingClaw we evaluated the grid dependency on the slump thickness after 600 s in a transect striking north–south. At the location of the thickest slump mass, the thickness obtained with Δ*x* = 185 m deviated only by 1.7% from that of Δ*x* = 93 m. The double the resolution of 370 m caused a corresponding deviation of 6.8%. Thus, we used a spatial resolution of 185 m. The slump was stable with a CFL number of 0.65.

We evaluated the spatial resolution of the input fluxes into GloBouss 120 s after failure by analysing the first wave amplitudes of the propagated waves along the same transect striking north–south. The wave height of the 926 m resolution differed only by 3.1% from that of 232 m. However, since it was feasible to use even 463 m in the modelling, we chose that. These sources were fed into GloBouss at various time intervals (see Table A2), whereas the resulting wave field 3000 s after failure was analysed. The amplitude of a resolution of 50 s deviated by 8% from a 30 s resolution. The wave amplitude of a finer resolution of 20 s deviated by 2% from the 30 s resolution. Thus, the flux sources were fed into GloBouss each 30 s.

We tested three spatial grid resolutions for the tsunami propagation model GloBouss 4100 s after failure. The first wave amplitude of the second finest resolution deviated 2% from the finest resolution, and the coarsest resolution deviated by 8% from the finest resolution. Since the finest resolution was feasible, we applied that one. The CFL number in GloBouss was chosen as 0.8.

### Kinematics

As we run our models in a depth-averaged regime, we divide the slump mass into vertical columns with length of one cell size. The height difference of the slump surface *H* in each column at two adjacent time steps serves as input for the vertical velocity calculation. With that velocity we calculate the vertical acceleration *a*_{z} in each column. It should be noted that resulting vertical velocities are half a time step behind the time of the surface heights of the next time loop, and the vertical accelerations are half a time step behind the velocities. The actual calculation of the vertical acceleration is a central derivative:
*t*^{(n = 1)} = 0.25Δ*t*, which means that the time interval for the calculations Δ*n* is not constant for *t* < Δ*t*.

Another kinematic quantity is the bed-parallel acceleration *a*_{||}, evaluated at the same time. In order to do so, we need to average the bed-parallel velocity *u*_{||} between two time steps and then evaluate the time derivative:
*n* = 1 applies here.

A third quantity is the angular momentum *m* is the mass of a vertical column, *r _{x}* and

*r*,

_{z}*v*and

_{x}*v*, respectively. The vertical velocity component corresponds to the one from the calculations above, but we approximate the horizontal velocity component

_{z}*v*with the bed-parallel velocity

_{x}*v*

_{||}, as the bed is nearly horizontal. Maximum bed slope angle is 5.25°. Equation (A4) shows the calculation for the total angular momentum, which is a sum of all angular momenta for each vertical column.

- © 2020 The Author(s). Published by The Geological Society of London

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/).