On 18 November 1929, an M_w 7.2 earthquake occurred beneath the upper Laurentian Fan, south of Newfoundland (Johnstone 1930; Hasegawa & Kanamori 1987; Bent 1995). The earthquake caused the largest observed historical landslide on Earth, displacing at least 100 km³ of sediment volume. This rapidly evolved into a turbidity current, as revealed by a series of successive breaks in telecommunication cables (Heezen & Ewing 1952; Heezen et al. 1954). A tsunami also occurred along the south coast of Newfoundland, causing a total of 28 fatalities along the Burin Peninsula (Doxsee 1948). The 1929 Grand Banks event is the only historical landslide-generated tsunami observed at transoceanic distances (Fig. 1), including locations along the east coast of the USA, Martinique, Bermuda, the Azores and mainland Portugal (Fine et al. 2005).

The first modern reanalysis of seismic signals suggested that the earthquake resulted from a major slump initiated just south of the Laurentian Channel, which rapidly spread laterally (Hasegawa & Kanamori 1987). They estimated a moment magnitude of M_w 7.2 (a seismic moment M₀ = 8 × 10¹⁹ N m). A more recent study by Bent (1995) suggested instead that the earthquake was due to a strike-slip double couple in a NW-striking plane, with a focal depth of c. 20 km. The location of the epicentre is shown in Figure 2. Figure 2 also shows a sketch of the initial slump area modified from the maps of Piper et al. (1985) and Hasegawa & Kanamori (1987). This area corresponds with the limit of instantaneous cable breaks (Fig. 3). It is thought that slumping and slope failures caused these breaks.

Piper et al. (1999) investigated the eastern part of the cable break area and found evidence for complex landsliding, including distributed slumps and translational landslides. The thickness of the sliding material was relatively shallow (c. 5–100 m thick, with an average thickness of 20 m; Piper et al. 1999; Mosher & Piper 2007). However, recently acquired seismic reflection data suggest that rotational slumping of a thick strata of unconsolidated sediment (c. 500 m) occurred on the St Pierre Slope, causing seafloor displacements (fault traces) up to 100 m in height. The slump occurred in the upper part of the slope at c. 1 km water depth, just 20 km east of the reported epicentre. The previously mapped surficial failures were thus a consequence of slumping of the thicker mass. The area of the fault is roughly A = 10 km × 20 km, which, combined with a low shear strength of μ = 10 GPa and a mean slip of about D = 50 m, is sufficient to account for the entire seismic moment M₀ = μAD = 8 × 10¹⁹ N m. Just south of the rotational fault from 1700 to 2300 m water depth, the new seismic reflection data show evidence of translational and retrogressive landslide scars. Previous interpretations by Piper et al. (1988), Piper et al. (1999) and Mosher & Piper (2007) indicate further additional failure along the western edge of the Grand Banks Valley drainage system and the Western and Eastern valleys of the Laurentian Fan. Our estimate of the entire failure volume is c. 135 km³ (Schulten et al. 2018). Further afield, successive cable breaks have left traces of a turbidity current moving with speeds of c. 15–30 m s⁻¹ (derived from Heezen et al. 1954), slowly decelerating as the masses moved away from the failure area and the slope angles decreased. According to Edgers & Karlsrud (1982), the mean landslide velocity may have been at least as high as 50 m s⁻¹ further up the slope.

A thorough review of the tsunami observations due to the Grand Banks landslide and earthquake are found in Ruffman (1997) and Fine et al. (2005) and some of their main findings are reviewed herein. In addition, recorded offshore water elevations, run-up heights and compilations of additional field investigations are found in the National Oceanic and Atmospheric Administration/National Center for Environmental Information global tsunami database (NCEI 2017). The 28 fatalities due to the tsunami occurred along the Burin Peninsula, but extensive coastal building damage was also reported in Cape La Hune on the south coast of Newfoundland. The maximum run-up heights along the Burin Peninsula were initially reported to be as high as 15 m (Johnstone 1930), later reconstructed to 13 m (Ruffman 1997). However, Fine et al. (2005) suggested that the run-up heights were significantly lower in places, ranging from 3 to 7.5 m. Back-analysis based on sediment transport models (Moore et al. 2007) suggest shoreline flow depth values of at least 7 m. This difference in the reports points to an appreciable uncertainty and a run-up range of c. 5–13 m seems plausible.

The wave along the Burin Peninsula first withdrew, followed by a breaking bore. Ruffman (1997) reported that the first wave arrival consisted of three pulses riding on top of a longer wave within 30 minutes. Further afield, the tsunami was noticeable in several places along Nova Scotia, with offshore water elevations at least as high as 0.6 m, but possibly as high as 2 m on the east coast (NCEI 2017). Wave gauge records in Halifax show that the period of the dominating wave was long, at least 30 minutes and even longer for the trailing wave system, with an amplitude of c. 0.6 m. Along the east coast of the USA, the observed offshore wave amplitudes range from 0.68 m in Atlantic City, New Jersey, to 0.15 m in Charleston, South Carolina. The wave caused destruction to boats in Bermuda and was clearly noticeable at transoceanic distances in Martinique, the Azores and in mainland Portugal (Fine et al. 2005), although records of tsunami heights are unfortunately not available from these locations.

Despite the importance of the Grand Banks tsunami, relatively few attempts to model the event exist. The most comprehensive approach was undertaken by Fine et al. (2005). They used a shallow water viscous fluid to model the landslide, coupling the landslide model to a shallow water tsunami model. The initial landslide consisted of an evenly distributed 5 m thick blanket at the start of the motion, covering the entire slump area shown in Figure 2. Visual inspection of their tsunami simulations suggests tsunami amplitudes just offshore Burin of 0.5–1 m, which are probably too small to explain the run-up observations. In the far field, simulated tsunami heights in Halifax were similarly 1 m, and c. 0.3 m offshore Atlantic City. Although their far-field simulation result seems to be roughly compatible with both the height and wave period observations for the US east coast locations, the directivity of the induced wave with similar amplitudes offshore Nova Scotia and Burin seems incompatible given the large contrast in reported run-up heights in these two locations. We also note that the wave traces simulated by Fine et al. (2005) in Burin start with a positive elevation, which conflicts with the eyewitness observations. Furthermore, Fine et al. (2005) did not report details of the computations (e.g. landslide speeds, grid resolutions) needed to dissect how the landslide emplacement process induced the tsunami. Much more effort is needed to understand how the Grand Banks landslide caused a large tsunami.

A second study by Trifunac et al. (2002) proposed that the tsunami was caused by a slump spreading rapidly in two directions. They used a Laplace transform approach to model the propagation of the wave. As their model apparently assumed a constant depth, wave directivity and refraction effects due to bathymetry were not considered, or at least not treated in sufficient detail. The induced simulated wavelengths were about three minutes, which may be reasonable to explain the short wave periods observed in Burin, but they are much shorter than the long-wave observations in Halifax.

This paper examines the landslide and the tsunami generation in more detail and discusses the source characteristics in light of both new geophysical observations of the source area and the tsunami and landslide observations. We provide more in-depth analysis of the simulated landslide kinematics and tsunami generation to shed light on the source processes necessary to induce the tsunami compared with, for instance, the most comprehensive previous study of Fine et al. (2005). To this end, the present work must be considered preliminary because there are still significant discrepancies between the observations and simulation results that need to be explained. On the other hand, the new modelling results may also rule out certain scenarios and show that the event was a complex one in which a series of local failures could be responsible for the wave propagating to the far field.

This paper is organized as follows: first, a brief review of the different numerical methods applied is provided. The study uses two alternative landslide models: (1) a prescribed slump model and (2) a translational viscoplastic landslide model. Both models are coupled to a linear dispersive wave model. In the results section, simulated slide motion as well as tsunami generation and propagation are shown while comparing the results with observations. We then discuss the results in a geological context and provide an hypothesis for the causes of tsunami generation.

Two independent source models are proposed for the landslide source. The first model is a pure rotational (rigid) slump, the motion of which is prescribed, whereas the second source model is based on a viscoplastic Hershel–Bulkley numerical model. In both cases, the simulated landslide seabed displacement is introduced as a time-dependent tsunami flux source.

The downdip landslide length is denoted with L, whereas the along-strike width is W. D is a measure of the thickness of the rotational landslide and ε is a shape parameter. The landslide moves along the x′ coordinate perpendicular to y′, according to a prescribed sinusoidal forward advancing travel distance function, travelling a distance 2R over a time period t_f = πt₀ (see Tappin et al. 2014 for details). Usually the slump source has a much smaller run-out distance R compared with the length L of the slide, which makes the source effectively mimic a rotational failure. For the tsunami simulations, we orient the slump source appropriately in the horizontal plane and interpolate it into the geographical coordinate system.

Here, γ is the total shear deformation. The equation system is first solved using a finite volume formulation for the conservative variables (LeVeque 2002) and then by a fractional step method.

In addition to the yield strength properties and hydrodynamic resistance properties, the bathymetry, the landslide density ρ_s and the Hershel–Bulkley parameter n are inputs to the model. The resulting outputs from the model are the instantaneous total landslide thickness (the sum of the plug layer and the shear layer) and the instantaneous velocities in longitudinal and latitudinal directions.

The GloBouss model (Pedersen & Løvholt 2008; Løvholt et al. 2008, 2010) was used to simulate wave generation due to the volumetric seabed displacement. The GloBouss model allows modelling of the tsunami propagation both by a linear shallow water model and by an optimized dispersive model. We use both versions of the model to illustrate the role of frequency dispersion on tsunami propagation. The seabed displacements were implemented as tsunami sources in terms of volumetric flux sources by replacing the ∂h/∂t terms in the continuity equation by a flux q subject to a low-pass filter (Glimsdal et al. 2013; Løvholt et al. 2015).

The slump source was placed close to a newly mapped rotational fault identified in the upper part of the St Pierre Slope (see Fig. 3; Schulten et al. 2018). Maximum vertical displacements of the slump within this fault zone range from ≤100 m in the northeastern part of the fault to ≤50 m in its western part. The thickness of the slump mass is c. 500 m. As shown in Figure 3, the main fault scarp strikes NE. The downdip length of rotational failure is c. 10 km and the along-strike length is c. 20 km. Estimates of the earthquake duration derived from seismic traces range from 15–40 (Bent 1995) to 50 s (Hasegawa & Kanamori 1987).

The fault orientations, displacements, lengths, widths and durations used in the slump model are compatible with these fault observations and, as shown in Figure 2, the modelled permanent change in slump displacement. The thickness D of the slump is set to 500 m, the shape parameter ε = 0.717 and the horizontal run-out distance of the slump centre of mass is 2R = 667 m, which gives maximum and minimum vertical displacements of ±100 m (Fig. 2). The slump volume is 17.5 km³, which is about three times larger than the 1998 Papua New Guinea slump source used by Tappin et al. (2008). The duration of slump motion was set to t_f = 60 s, i.e. similar to the estimates of Hasegawa & Kanamori (1987) for the earthquake duration. The downdip length and along-strike width were both set to 10 km. Initial simulations used a width of 20 km, but this dimension caused a large tsunami with large amplitudes over a greater coastal stretch than just the Burin Peninsula, and was hence incompatible with the tsunami observations. The slump orientation was roughly aligned with the observed fault. In the slump model, the x′ direction pointed northwards and the source was rotated 200° counterclockwise to roughly comply with the orientation and direction of mass movement of the slump.

Tsunami simulations were carried out using the slump source parameters, tsunami propagation models and source filters described above on a bathymetric surface refined from ETOPO 1 (www.ngdc.noaa.gov/mgg/global/) with a grid resolution of ¼′ × ¼′. The fine grid resolution was necessary as the short duration and horizontal length scales of the slump resulted in short tsunami wavelengths. Using a coarse grid (e.g. 1′ × 1′) the waves would not propagate over the relatively shallow shelf because they would not be sufficiently resolved. A Courant number of 0.8 was used. The waves were strongly dispersive, hence the simulations were carried out using an optimized dispersive version of the GloBouss model. The landslide fluxes were fed into the tsunami model every 30 s.

The maximum surface elevation due to the slump source is shown in Figure 4. The surface elevation decays rapidly southwards and azimuthally from elevations of several tens of metres near the source to a few meters 100 km away from the landslide. The short lateral extension of the source favours strong radial spreading. Northwards, the bathymetry guides the wave towards the Burin Peninsula, other parts of southwestern Newfoundland and St Pierre and Miquelon. Large near-coastal sea surface elevations up to 10 m are found in the tsunami simulations. As elaborated in the folllowing, the orientation of the slump motion also influences the wave radiation pattern. The short dispersive waves in combination with the shallow bathymetry give rise to a complex refraction pattern.

The sensitivity of the slump orientation was investigated, carrying out simulations with horizontal fault orientation angles of 160° (a slight westwards orientation) and 180° (a southwards orientation). The sea surface elevation time series resulting from the simulations offshore the Burin Penninsula (Lawn and Burin) and offshore western Nova Scotia (Glace Bay) are shown in Figure 5. These time series gauges show that the tsunami wave period is two to three minutes. With respect to eyewitness observations reporting three short waves appearing within a time frame of 30 minutes, the wave periods are rather short, but the time series clearly also show fluctuations implying that not all the peaks in the wave train were equally strong. The maximum offshore wave amplitudes in Lawn and Burin for the 200° source range from c. 5 to 6 m and may be further amplified due to shoaling and coastal effects. These coastal effects were not sufficiently resolved by the offshore model and more exact estimates of the run-up emerge from coupling the offshore wave propagation to a coastal inundation model (e.g. Løvholt et al. 2010). On the other hand, the time series points are placed close to the shoreline and therefore a significant portion of the shoaling should already be captured. Offshore Glace Bay, Nova Scotia, the surface elevations due to the 200° source do not exceed 1 m. The two other source orientations, and particularly the 160° source, give a more equal distribution of maximum wave amplitudes, with smaller amplitudes offshore Burin and larger amplitudes offshore Nova Scotia. Hence aligning the slump source along the fault does not only comply with the new geophysical data, but it is also favourable with respect to the relative differences in wave observations.

Sensitivity studies with respect to the slump position and duration (and maximum speed and acceleration) were also conducted, although the results from these studies are not shown. Simulations with other fault positions (position shifted 9′ eastwards) did not cause systematically higher or lower waves than those shown in Figures 4 and 5. However, the simulations with a longer duration of slump motion (t_f = 120 s) resulted in roughly halved wave amplitudes (compared with t_f = 60 s) and waves that were smaller than the field observations. Simulations with shorter slump durations (e.g. t_f = 30 s) also produced excessively high tsunami amplitudes.

Simulations coupling the translational landslide model to the tsunami propagation model were conducted. The landslide simulation study area roughly covers the area depicted in the lower panel of Figure 2, whereas the tsunami simulations were carried out over the domain shown in Figure 1. Both the linear shallow water model and the dispersive version of GloBouss were used in the tsunami simulations. Grid resolutions were ¼′ × ¼′ for the landslide model and 1′ × 1′ for the regional tsunami simulations.

The initial landslide geometry is shown in Figure 2. Compared with the slump source, the translational slide masses are thinner, but cover a much larger area. In the initial landslide configuration, the landslide-prone masses are placed between water depths of c. 500–2500 m, with an initial maximum thickness of 40 m in the proximity of the continental shelf, but with linear thinning as the depth increases (giving a mean thickness of 20 m). A series of initial model runs with a volume of c. 100 km³ and small, yet realistic, flow resistance parameters provided tsunami heights that were almost an order of magnitude smaller than the reported far-field wave observations.

To improve the consistency of the model results with the actual wave records, we first had to use a larger slide volume and then place a larger fraction of the slide masses in shallower water (source parameters are given in Table 1). The new slide volume used was 260 km³, exceeding our best estimate from field investigations, which ranges from 140 to 180 km³. Although the applied flow resistance parameters are low and the landslide volume is high compared with the field observations, the resulting landslide kinematics serve to demonstrate the critical dimensions for appropriate tsunami generation in terms of extent, volume, speed and water depth. We will attempt to shed light on the wave generation process, and the most important factors determining wave generation, by visualizing the effects of different tsunami source strength parameter values.

Figure 6 shows the simulated maximum surface elevation due to the translational landslide. The maximum elevations are not as high as for the slump source shown in Figure 4, but this source produces larger wave amplitudes over a much wider area. The complex bathymetry creates a complicated refraction pattern. It is difficult to see a dominant wave directivity from Figure 6.

The time series of the surface elevations close to the Burin Peninsula generated by the translational slide simulation show water elevations of 2.5–3.5 m at several coastal locations (Fig. 7). The wave first withdraws and then rises. The wave arrival consists of a series of short waves overriding a longer undulation and these waves inherit longer wave periods than the waves induced by the slump source. Both of these wave characteristics are consistent with eyewitness observations. In Halifax, the wave starts to rise and then causes a trailing wave system with wave periods of c. 20–30 minutes. Importantly, the simulation results for Halifax are consistent with the main characteristics of the observed Halifax wave trace (see Fine et al. 2005). First, the simulation more or less replicates the shape and polarity of the first wave pulse; second, it produces wave periods that are at least partly consistent with the observed wave trace; and third, the simulations produce the correct positive wave amplitude of 0.6 m (corrected for tide elevation by Fine et al. 2005). For the other location offshore Glace Bay, Nova Scotia, we obtain wave elevations >2 m, which is slightly higher than expected based on eyewitness descriptions from various locations along the east coast of Nova Scotia (see NCEI 2017 for descriptions). For the three coastal locations offshore the eastern US coastline (New Jersey, Maryland and South Carolina), our simulations are close to the surface elevations from reported wave records (see Fig. 1 for the observed heights). As noted earlier, some additional amplification at coastal points is expected (see Løvholt et al. 2017 for analysis of wave amplification due to the Storegga Slide tsunami). In Bermuda, wave elevations up to c. 0.7 m are simulated. Bermuda is the only one of the investigated gauge points that are strongly influenced by frequency dispersion.

Observed and simulated tsunami travel times (Table 2) show satisfactory agreement between the simulated and reported wave arrivals because the estimated arrival times are probably a mix of observations of the first and peak arrivals. In addition, the actual arrival times are sensitive to the bathymetry, whereas the simulation uses a coarser grid representation of this bathymetry. Importantly, this comparison substantiates the view that wave generation must have occurred almost immediately after the earthquake and not as a result of the far-field turbidity currents recorded south of the landslide release area several hours after the earthquake (Heezen & Ewing 1952).

Figure 8 shows the initial stages of the simulated landslide and the contemporaneous tsunami generation. Snapshots 5, 10, 20 and 30 minutes after the release of the landslide are shown. After five minutes, a chaotic picture with several slope failures in different directions is seen. The slope failures initiate in steep sections of the bathymetry, focusing towards the canyons. This focusing leads to increased thickness of the landslide. Despite the chaotic pattern of slope failures deflecting the flowing masses in different directions, the initial wave has a distinct shape with a pronounced elevation forming downslope and a depression upslope. After ten minutes, the major part of the mass movements takes place in the canyons. Until this point, the wave build-up has increased, causing most of the wave energy to propagate towards the SE. In addition, a secondary wave moving SW forms in the eastern part of the slope area. After 20 minutes, larger structures are formed from different branches of the landslide. At this time, the majority of tsunami generation has taken place. After 30 minutes, a clear landslide front forms as the landslide reaches a gentler part of the slope. The significant part of wave generation has passed at this stage, evident from visual inspection of the simulated wave evolution and the match between simulated and observed wave arrival times.

Within about one hour the landslide has reached the locations of the first cable breaks downstream (results not shown). The simulated landslide run-out does not reach the locations of the cable breaks that occurred after about three hours. It is widely understood that the sequential cable breaks were a consequence of a turbidity current rather than of the initial landslide and consequent debris flow (Heezen & Ewing 1952). On the other hand, such turbidity currents can be eroded from the top of the initial dense flow, in particular for flows involving high sand contents (Ilstad et al. 2004; Elverhøi et al. 2005; Breien et al. 2010). In fact, rapid flow transformation of the initial landslide into debris flows and then channelized turbidity currents is indicated (Piper et al. 1999; Mosher & Piper 2007). In this model, a landslide motion over the first c. 10–20 minutes would be sufficient to generate the observed tsunami, whereas the later stages of landslide motion influence tsunami generation to a smaller extent.

Here, u and v are the velocity components of the landslide, g is the gravitational acceleration, and h is the water depth. The second quantity is the product of the instantaneous flux q multiplied by the Froude number (q × Fr). As noted earlier, the flux q represents the landslide source strength that is input into the tsunami model at a given time step. The Froude number measures the speed of the landslide relative to the wave and hence the criticality of the wave generation. A unitary Froude number means that the landslide and the tsunami move with the same speed, which renders an efficient wave generator. In turn, smaller Froude numbers imply less efficient wave generation (e.g. Ward 2001; Harbitz et al. 2006; Løvholt et al. 2015). Consequently, the product of the flux and the Froude number is a measure of the wave generation efficiency. In this model, the local Froude number extends up to c. 0.7, which is a high value for deep water submarine landslides. We see that the largest Froude numbers occur at 10 and 20 minutes and in the shallowest part of the slope between 500 and 2000 m water depth. The Froude number is significantly reduced after 30 minutes. The wave generation potential of the Grand Banks landslide becomes even more evident when comparing the q × Fr product with the wave propagation snapshots. This quantity rapidly increases during the early stages of motion, reaching a maximum at about ten minutes, and is clearly reduced after 30 minutes.

Figure 10 shows the product q × a/c₀, which can also be taken as the temporal derivative of q × Fr. In studies of other landslides, Løvholt et al. (2005) found that the product of initial acceleration and landslide volume governed tsunami generation for a block landslide; Haugen et al. (2005) suggested that for the source strength, the acceleration should be scaled by the local wave celerity. Compared with q × Fr, q × a/c₀ seems to attain slightly larger values in the early stages of motion and reduces more rapidly. Because our inspection of the wave propagation fields suggests that the main wave generation has taken place within 20 minutes, q × a/c₀ seems to be a good indicator of the landslide source strength. It should be noted that the scaled initial acceleration is considered to be more important than the scaled acceleration at later stages of motion because rapid initial acceleration prevents an early cut-off of the front and rear parts of the initial waves, which allows the waves to build up amplitudes (Haugen et al. 2005). In any case, the high values of q × Fr and q × a/c₀ during the first 10–20 minutes both correspond with the time when the tsunami generation is strongest. The link between tsunami generation and Fr shows that rapid movement and a high acceleration of mass in the shallowest region of the slope is a necessary condition for the effective generation of a tsunami by of the Grand Banks landslide. A landslide occurring at larger water depths would not be sufficiently tsunamigenic.

We used two different models, the first a pure slump model and the second a translational landslide model, to simulate the generation and propagation of a tsunami by the 1929 Grand Banks landslide. The slump model complies with both new geophysical data showing a large rotational failure and the duration of the fault motion from seismic records. The slump-induced tsunami simulations explain the near-field observations of the waves in the Burin Peninsula, Newfoundland, but not the long period far-field waves observed in Halifax, Nova Scotia. The simulations are not sufficiently extensive to explain the wide range of instantaneous cable breaks that occurred in the Laurentian Fan. The translational landslide source covers the area of instantaneous cable breaks. The simulated waves provide a good overall match with wave observations in both the near and far fields, with wave amplitudes possibly too high towards the east coast of Nova Scotia, but possibly lower than observations close to the Burin Peninsula. Placing the landslide masses in a shallow depth of water (c. 500–1500 m) was necessary to produce sufficiently large waves in the far field. Piper et al. (1988) suggested that the area outlined by the initial cable breaks (Fig. 3b) represents the area of instantaneous landsliding. Piper et al. (1999) show evidence of landsliding at the heads of the Laurentian Fan valley systems (e.g. the Eastern Valley), which is c. 500 m long. The observation of long waves is only compatible with a landslide event of large areal extent.

The simulation results herein argue that a combination of a slump and widespread translational landsliding could have caused the Grand Banks tsunami. The slump provides the most likely source for the large run-up of the tsunami observed along the Burin Peninsula, whereas widespread translational landsliding is responsible for the far-field tsunami. The Grand Banks landslide was undoubtedly a complex event with interactions between many different slope failures causing a large tsunami. We carried out preliminary simulations that provide hints of some possibilities and constraints. Many other landslide volume configurations and dynamic parameters need to be tested. Refined simulations with a higher coastal resolution can be carried out near the shore to obtain more accurate results (e.g. Løvholt et al. 2008; Tehranirad et al. 2015). We note that the simulations are associated with large uncertainties, both with respect to the initial failure and simplification in the slide representation, such as the depth averaging. We plan to unite the simulations with detailed geophysical observations in future analyses to give a more coherent description of the Grand Banks landslide and tsunami.

The work was funded by the Research Council of Norway project Tsunamis Induced by Large Landslides (NFR 231252/F20). The funding agency was Norges Forskningsråd (ID 501100005416). The principal award recipient was Finn Løvholt.