Rheological considerations for the modelling of submarine sliding at Rockall Bank, NE Atlantic Ocean

The post-failure stage of submarine landslides can be adequately described using the fluid mechanics principles and is often considered as the most significant for tsunami generation.^1,2 Above all, the sediment concentration and water content of the landslide play an important role on the rheological behavior. For example, landslides with water content greater than 50% behave as Newtonian fluids (classified as streamflows), whereas landslides with less water content exhibit a visco-plastic rheological behavior (classified as slurry and granular flows).³ 

Depending on the selection of the input parameters for the calculation of the shear retarding stress, the flow properties may vary corresponding to different fluid types. Three main types of rheological models have been used in the past to describe the dynamics of submarine landslides (considered as non-Newtonian fluids) at the post-failure stage:^1,4 (1) viscous models, (2) frictional models, and (3) visco-plastic models. De Blasio et al.⁵ referred to the latter two as the main flow models for the simulation of submarine motion.

The appropriate rheological regimes that can be used in the simulations may vary between different environmental settings depending on various factors. Among those are the sediment composition, the particle size distribution and shape, and the pore pressure.^6,7 For example, three types of fluids have been commonly used to represent the rheology of mudflows or muddy debris flows (silt-rich and clay-rich flows): a Bingham fluid, a Herschel-Bulkley fluid, and a bilinear fluid.¹ On the other hand, terrestrial events and sand-rich flows exhibit a granular behavior where there is a lot of grain-to-grain interaction. These types of flows are better described by frictional flow models.^5,7,8

Following the description of the rheological regimes as found in Ref. 2, the Bingham rheology yields

where T is the magnitude of the retarding stress of the flow, T₀ is the yield strength, representing a critical threshold beyond which the body begins to flow, μ is the dynamic viscosity, and $γ̇$ is the shear rate (or also strain rate). For zero viscosity, the second term is neglected and the flow motion becomes subject to a constant retarding stress.⁹ 

The Herschel-Bulkley rheology yields

where K is a consistency factor and n is the flow index. For n < 1, the state of mixture describes a shear-thinning fluid (or else a pseudoplastic fluid), a fact which indicates that the viscosity decreases as the shear rate increases. On the other hand, for a shear-thickening fluid (or else dilatant fluid) where n > 1, the viscosity behaves in the opposite way from that of a pseudoplastic fluid; this phenomenon is rather uncommon and is mostly observed in sand-water mixtures.¹⁰ In the case where n = 1, the formulation reduces to that of a Bingham fluid. Finally, a bilinear fluid which can be described as a generalization of a Bingham fluid⁶ is given by

where γ₀ = T₀/(μ_h − μ_l). The existence of two different dynamic viscosities μ_l and μ_h highlights that when the shear rate is high the flow behaves as a Bingham fluid with low viscosity (μ_l), but at low shear rates, the flow behaves as a Newtonian fluid with high viscosity (μ_h).⁶ 

The above rheological models have been commonly used for the numerical modelling of cohesive submarine flows.^5,6,11 Imran et al.⁶ have tested the suitability of the three models and observed that each model influenced differently the run-out length of the flow, for example, a bilinear model resulted in a larger run-out distance than a Herschel-Bulkley model. The authors concluded that the bilinear model gives a better representation of the stress-strain relationship, particularly at low shear rates.⁶ The sediment composition can also play an important role for the choice of the rheological model. The Herschel-Bulkley and especially the Bingham rheology are very common for the modelling of muddy debris flows.^5,7,12

On the other hand, for the description of granular flows, the Mohr-Coulomb failure criterion is often used. The rheology of the failing mass (see also Ref. 13) can be expressed as

where c is the cohesion of the material (for a purely frictional material c = 0), σ is the stress acting normal to the surface, ι is the pore water pressure, and ϕ is the Coulomb friction angle. The term (σ − ι) represents the effective normal stress. The flow resistance to motion is a result of the frictional forces dominant between the particles of the mass and the basal layer of the flow. A linear relationship exists between the shear and the effective stress. Three types of shear stresses can be considered the principal mechanisms driving material down-slope for submarine landslides: gravity loading, seismically-induced stress, and storm wave-induced stresses.¹⁴ Gravity driven flow behavior is often incorporated in the landslide modelling.

Earlier studies on the mechanics of submarine landslides indicated that both viscous and plastic terms become important on the simulation of complex flow motion.¹⁴ Norem, Locat, and Schieldrop¹⁵ incorporated both visco- and plastic-properties in the rheological regime. The retarding stress could be thus computed (see also Refs. 2 and 14) by

where r_ι is the pore pressure ratio equal to ι/γh, ι is the pore pressure, and γ = ρg is the specific weight of the sediments with h thickness. The approach in Ref. 15 stresses the importance of incorporating a drag term to account for the forces exerted by the seawater on the frontal and upper body of the flowing mass. The incorporated water drag can reduce significantly the maximum velocity of the landslide.

To account for viscous drag forces, Hutter and Nohguchi¹⁶ made use of a viscous Voellmy-type resistive stress term, proportional to the square of the landslide velocity. More recently, Løvholt et al.¹⁷ have combined a turbulent term with the granular model introduced by Savage and Hutter¹³ to incorporate the water drag in their submarine landslide simulations. The turbulent term results from the empirical formula first described by Voellmy for the rheology of snow avalanches.¹⁸ The formula can be written as

where T is the stress vector, u is the flow velocity, and ξ is a non-dimensional coefficient used to represent the effect of turbulence and/or grain collisions.¹⁹ 

As submarine flows are almost impossible to observe directly, the rheological properties of the flow can be generated by geotechnical measurements in the laboratory. Laboratory experiments are also used to understand better the dynamics and determine experimentally some values; they can also be used in conjunction with numerical experiments for validation purposes.^{11,13,20–22} However, the values of the parameters may not be always easily defined due to the lack of data/facilities, etc. In the absence of adequate input data, a common practice in the numerical modelling of submarine slides is to try to match the depositional lobes resulting from the simulations, with the landslide deposits observed in the field. A drawback in this approach is that the simplifications of the physics in the models cannot capture the complex physics of the landslides.

For example, a constant retarding stress can be often assumed for the whole duration of the flow whereas this is rarely the case in reality. Mixing of the slide material with water increases its pore water pressure and can have a significant reduction in the shear strength of the flow.⁵ Measurements of the yield strength on the sediments of the slope and on the deposits can thus differ substantially.²³ In addition, mechanisms such as hydroplaning and remoulding that enhance the capability of the flow to travel large distances are often not accounted in the modelling.

Lastras et al.²⁴ have found that for the simulations of the Big’95 debris flow, the best fit solution had a very low yield strength of 800 Pa. The authors have also noted that the internal friction had also a dramatic influence on the run-out length of the flow and best fit values were as low as ∼1°. They concluded that the absence of hydroplaning in the modelling could explain the very low yield strengths required. Simulations of hydroplaning and non-hydroplaning debris flows in Ref. 7 show that under a non-hydroplaning regime the yield strengths have to be significantly reduced for the flow to reach the same run-out distance.

The sediment composition of the slide material has an effect on the cohesion of the flow. Sediments with high clay content suggest more cohesive landslides and can thus be described by higher yield strengths (such is the case for sliding in the Bear Island Fan).¹¹ On the other hand, mixtures of clay and sand gravel appear to be less cohesive (like the 1929 Grand Banks slide¹¹). Increase in the yield strength of the flow is linked with decrease in the observed run-out distance. De Blasio et al.⁵ have observed that an increase in the yield strength from 6 to 14 kPa in the simulations resulted in a decrease in the run-out distance of the flow by a factor of 2. The value of the yield strength can thus be of great importance for the modelling of submarine landslide deposition, and it may vary between different numerical models.

In the present work, we model the episodes of collapse at the Rockall Bank Slide Complex (RBSC), NE Atlantic Ocean under a rheological framework. We note that the terms landslide and slide in this case refer to general slope collapse in the region without implying any particular type of motion. For the landslide simulations, the one-fluid version of the two-dimensional code VolcFlow has been employed.⁹ VolcFlow is a finite-difference Eulerian code that solves the general shallow water equations of mass and momentum conservation. The algorithm can account for a simplified representation of varying rheological regimes. It has been efficiently used to study many cases of subaerial and submarine sliding.^19,25,26

A frictional model, based on the Mohr-Coulomb failure criterion, was used to study the flow motion. A simplified Bingham rheology is also implemented. The results of our simulations show that the Bingham rheology represents successfully the landslide deposition for this case, whereas the basal friction critically affects the run-out length of the flow and thus has to be significantly small. To study separate failure episodes, we make use of a model that combines the Bingham rheology with frictional terms for the rheology of the landslide. In all the scenarios, a velocity dependent term is incorporated in the models to consider the effect of hydrodynamic drag in the motion. We make use of a Cartesian coordinate system, where the x- and y- axes represent the EW and NS horizontal directions, respectively, and the z- axis represents the vertical direction. The retarding stress $T=Tx,Ty$ is given by

where du/dh is the shear rate, r is the slope curvature, and ϕ_bed is the basal friction angle.

As it can be observed, the first three terms of the equation resemble the model introduced by Norem et al. to describe the mobility of subaqueous flows.¹⁵ The last term of the equation is making use of the coefficient of turbulence; in VolcFlow, the term is non-dimensional.¹⁹ In this paper, we first present an overview of the Rockall Bank and multiphase failure in the area in Sec. I B. The methodology of the work (Sec. II A) and the results of the varying landslide scenarios (Sec. II B) are given in Sec. II. A short discussion on the tsunamigenic consequences of the landslide scenarios is also given in this section (Sec. II C). Finally, the conclusions of this work are presented (Sec. III).

The Rockall Trough is a steep-sided sea-floor depression approximately 1000 km long and 200–300 km wide. It forms the largest of the three basins located northwest of the UK and Ireland in the Atlantic Ocean, the other two being the Hatton Basin and the Hatton Continental Margin²⁷ (Fig. 1). It is also the closest of the three in proximity to the Irish coast.

The trough is bounded on the east by the Irish and Scottish continental shelves; Rockall Bank, an almost flat-topped, underwater plateau lies to the west (Fig. 1). The bathymetry of the Rockall Trough ranges between 300 and 400 m on its east and west margins, abruptly going down to 3000 m at the north part of the basin, deepening gently to 4000 m toward the south, at the Porcupine Abyssal Plain (Fig. 1).

Two orthogonal wide-angle seismic profiles, recorded during the Rockall And Porcupine Irish Deep Seismic (RAPIDS) project (1988-1990) together with previously attained gravity data revealed that the Rockall Trough formed as a result of a multiphase rift evolution that took place in the Triassic, Late Jurassic, and Early Cretaceous.²⁷ In more detail, the data revealed the existence of continental crust lying beneath the three major sedimentary basins: the Rockall Trough, Hatton Basin, and Hatton continental margin.

Three major depositional areas are distinguished on the seafloor of the Rockall Trough: the north part is occupied by the Donegal-Barra Fan and the deposits of the Rockall Bank Slide Complex (RBSC), also referred to as the Rockall Bank mass flow (Fig. 1). The areas lie offshore the Irish and Scottish continental shelves. The third one, the contouritic Feni Drift, extends on the west margin of the Rockall Trough and up to the deepest part of the basin in the south^29,30 (Fig. 1). Among the three, the RBSC constitutes the largest region of submarine slope failure in the Irish Atlantic margin, affecting a total area of 18 000 km².²⁹ 

Roberts was the first to describe the submarine failure scarps that extend widely on the eastern slopes of the Rockall Bank.³¹ Several bathymetric surveys have been carried out in the area ever since. The Geological LOng Range Inclined Asdic (GLORIA) side-scan sonar survey (the Atlantic Irish Regional Survey—AIRS96) took place in 1996 and the Towed Ocean Bottom Instrument (TOBI) survey of a higher 5-10 m spatial resolution in 1998.^32,33 Between 2000 and 2001, R. V. Bligh collected multibeam bathymetric data of most of the Irish offshore area as part of the Irish National Seabed Survey (INSS) program. The multibeam bathymetry and side-scan sonar data acquired during the surveys have provided a full coverage of the slope failure features.

Two distinct areas of failure escarpments (upper and lower slope regions) have been identified on the eastern margin of the Rockall Bank.²⁸ Due to its diverse scarp morphology, the upper slope was further subdivided into three distinct parts: the north, the south, and the central upper slope.²⁸ The two regions are separated by an elongate deep erosional moat near the base-of-slope, which traverses the area for 120 km and in parts can be as deep as 150 m.^28,29 The gradients of the slope are steep varying between 5° and 10° and in places can reach up to 15°-20°.²⁸ The scarps vary on sizes and heights. Despite some isolated retrogressive failure episodes, an overall non-retrogressive failure mechanism has been proposed based on observations of the scarp morphology and geometry.²⁹ 

Downslope of the scarp complex is the depositional area occupied by sedimentary lobes (Fig. 2). The deposits of the RBSC show a run-out length of 120 km and a depositional width of 150 km^28,29 [Fig. 2(c)]. The width to length aspect ratio may be characterised as irregular when compared to other submarine landslide deposits found on the western and eastern margins of the Atlantic Ocean.^29,34 The occurrence of turbidites in piston cores, collected from the deepest part of the Rockall Trough basin, supports that the flow might have transformed into a more-dilute suspended mixture during motion that had travelled as turbidity current.³⁵ Sedimentary core data suggest a mixture of clayey muds, silt, and sandy layers for the composition of the flow.³⁵ 

In the past, there has been a high degree of uncertainty regarding the triggering mechanism of the failure and whether sliding occurred as a single failure episode or multiple successive or separate episodes of failure.²⁸ One hypothesis for the triggering mechanism suggests that sedimentation on top of abrupt pre-existing fault scarp topography at the upper part of the slope enhanced by strong bottom current activity removed support from the base-of-slope, leading to slope instability and subsequently failure.^28,29 Elliott et al. concluded that the contourite drift dynamics in the basin acted as a substantial long term factor promoting slope instability.²⁹ The sedimentation pattern has caused overloading of the upper part of the slope, whereas the contourite drift systems contributed to the erosion of the lower part of the slope and led to the formation of an elongate moat. Georgiopoulou et al. found a direct relationship between the RBSC seafloor scarps and the basement morphology and suggested that focused fluid flow along basement-bounding faults and/or differential compaction across the scarps contributed to the sliding.²⁸ 

A manual reconstruction of the preslide morphology estimated that the volume of the missing sediments ranges between 265 km³ and 765 km³ (conservative and generous approach, respectively, in Ref. 28). Sedimentary and seismic data, acquired from the depositional lobes of the RBSC during the RV Celtic Explorer cruise CE11011, support the hypothesis that the RBSC morphology is the result of multiple phases of slope failure.³⁵ At least three major phases of collapse can be distinguished based on the data.^28,35 In more detail, the seismic profiles along the lateral extents of the lobes have revealed the limits of three depositional horizons separated by long periods of slope stability³⁵ (Fig. 2). The two first episodes of collapse (here referred as landslide A and B) have more likely initiated from the south upper slope (SUS) region of the margin, whereas the most recent episode (landslide C), the lobes of which are still visible on the seabed of the Rockall Trough, has resulted from collapse on the north upper in combination with failure from the lower slope region.³⁵ 

As landslide A forms potentially the oldest failure episode and is buried deeply, its age is hard to estimate; it is suggested that it must be a few Ma old.³⁵ An estimation of the average sedimentation rate on the slope scarps suggests that landslide B may have happened at ca. 200 ka.³⁵ Interpretation of radiocarbon dating shows that the most recent event took place 21 ka, which is 3000 yr post Last Glacial Maximum (LGM), when the British Irish Ice Sheet (BIIS) was still at maximum extent and was only just starting to destabilise.³⁵ 

The volumes of the three phases (landslides A, B, and C) were estimated based on information from the sediment piston cores and the 2D seismic profiles of the landslide deposits at the complex.³⁵ The volume of landslide A was estimated at approximately 200 km³, the deposits cover an area of 7500 km² [Fig. 2(a)].³⁵ The second phase (landslide B) was smaller in size occupying an area of 4500 km² with a material volume of ca. 125 km³ [Fig. 2(b)].³⁵ The data regarding the 3rd failure episode (landslide C) yield a volume of ca. 400 km³, occupying an area of 6600 km² [Fig. 2(c)].³⁵ 

Following these volume estimations, the total volume of the deposits becomes 725 km³. This comes in close agreement with the estimated maximum volume of the missing sediments from the slope (765 km³), i.e., the generous approach for the reconstruction of the preslide morphology.²⁸ It is noted that the two older episodes (landslides A and B) involve higher levels of uncertainty than the most recent one (landslide C). The estimations of the depositional extents of landslides A and B are bounded by the limited number of the seismic profiles acquired during the surveys. The event for which the most information exists is the 3rd episode of collapse, the extent of which is visible on the seafloor of the trough. As a consequence, the landslide modelling of the events draws information from the latest episode of collapse to study the landslide deposition.

For the landslide modelling, we follow closely the methodological approach described in Ref. 36. The bathymetry data were retrieved from the EMODnet Bathymetry portal. The data were transformed to the Cartesian coordinate system, smoothed, and interpolated to form a uniform grid. A Digital Elevation Model (DEM) with dimensions of 200 km × 250 km and spatial resolution Δx = 320 m, Δy = 500 m was employed for the landslide modelling. VolcFlow computed the run-out length of the landslide. To establish the deposition, the landslide motion was initially modeled with a time step that varied in accordance with the stability criteria in the code. A stopping criterion for a velocity lower than 4 m/s was used, as no significant depositional re-arrangement was observed beyond that point.

For each scenario, a two-dimensional Gaussian function is used to represent the initial body of landslide. The Gaussian shape forms a simplified yet efficient and standardised approach to approximate the landslide body in numerical modelling. Landslides of a Gaussian shape exhibit a smoother mass distribution and integrate well with the surrounding topography resulting in a better representation of the pre-slide morphology.^21,37–39 In the numerical simulations of submarine sliding in the RBSC, a landslide of a Gaussian shape was computed with the mass spreading along the x- and y- horizontal directions (Fig. 3). The shape of the mass is given as a function of the landslide thickness h and the standard deviation of the landslide, σ_x and σ_y.

Note that Sarri, Guillas, and Dias⁴⁰ performed a sensitivity analysis of the Gaussian shape impact on the tsunami wave created. In that study, there was no significant sensitivity to the spread ratio (i.e., the ratio of the characteristic length over the characteristic width for a constant volume). However, this may be due to the smooth aspect of a Gaussian shape: other irregular (and time-varying) shapes may have a different influence on the resulting wave. However, the investigation of the influence of a high dimensional shape as input is complex and was only recently addressed⁴¹ in the context of the influence of bathymetry onto tsunami wave propagation.

The reconstruction of the most recent episode (landslide C) was performed as more information exists for its depositional lobes.³⁵ Taking under consideration the best-fit values of the rheological parameters for landslide C, the numerical simulations for landslides A and B were performed. This modelling approach follows the hypothesis that the sediment composition in the area is fairly similar. As introduced in Sec. I B, the slope failure that has resulted in the depositional lobes of the last episode possibly initiated from the North Upper And Lower Slope (NUS-LS) regions. The scarp heights on the slope region are measured to range between 200 and 300 m at the north upper and south upper slopes. The lower slope is dominated by smaller scarp heights that range between 60 and 120 m.

These scarp heights were used for the initial thickness of the landslide. A landslide with maximum thickness $hmaxNUS=300 m$ and a volume of V_NUS = 264 km³ was assumed for the north upper slope. For the lower slope region, a landslide of $hmaxLS=120 m$ and volume V_NUS = 136 km³ was considered (Fig. 3). We took under consideration the maximum scarp heights of the landslides first to account for a worst case scenario and second to avoid excessively large spatial extents of the landslide on the slope. Finally, the total volume of the event was V₃ = 400 km³ matching the estimated volume of landslide C. For the older episodes, we consider two landslides initiating from the south upper slope region (landslides A and B). Landslide A corresponds to the oldest recorded failure episode in the region with an estimated volume of ca. 200 km³.³⁵ Another collapse initiating from the same slope area with a volume of approximately 125 km³ (landslide B) happened at a much later stage.³⁵ 

Table I presents landslide scenarios examined with the three rheological models. Seven landslide scenarios are presented here, of which scenarios 1–5 simulate landslide C and scenarios 6 and 7 simulate landslides A and B, respectively (Table I). Scenarios 1 and 2 make use of a purely frictional rheology similar to Eq. (4), scenarios 3 and 4 represent a Bingham rheology with frictional terms [described by Eq. (7)], and scenario 5 studies the effects of the Bingham rheology on the landslide motion [Eq. (1)]. The main rheological parameters incorporated in the numerical simulations are the yield strength of the landslide, T₀, the basal friction angle, ϕ_bed, the dynamic viscosity, μ, and the coefficient of turbulence, ξ, used to represent the hydrodynamic drag. The selection of the appropriate values of the rheological parameters is based on the literature and the results of the Bayesian calibration in Ref. 36.

The one-fluid version of VolcFlow assumes that no mixing is taking place between the water and the landslide. To partially account for mixing, a reduced density is used for the landslide (ρ = ρ_ls − ρ_$w$, where ρ_ls denotes the landslide density and ρ_$w$ denotes the water density, respectively). In numerical modelling of submerged motion, such an approach has been used for many cases.^20,26,42 Here, we consider a fixed density of ρ = 1, 200 kg m⁻³.

Landslide scenarios 1 and 2 model a frictional flow motion with the incorporation of a velocity dependent term [Table I and Figs. 4(a) and 4(b)]. It is observed that the basal friction has to be very low to produce a run-out length [indicated with a black dashed line in Fig. 4(a)] that matches with the one of the field deposits (orange dashed line). In more detail, a flow that has ϕ_bed > 0.5° results in deposits that do not move significantly but rather accumulate close to the slope region. They exhibit a more lateral than longitudinal spreading [Fig. 4(a)]. The basal friction has to be significantly reduced down to ϕ_bed = 0.2° to produce a run-out length that matches the observations [scenario 2, Fig. 4(b)]. However, in this case, the landslide deposits exceed the lateral limits of the depositional lobe and also exhibit large spreading close to the slope region [Fig. 4(b)].

It is noted that the dimensional coefficient ξ in both scenarios is equal to 0.045 (Table I). The term has a direct impact on the maximum velocity of the landslide which rises up to 60 m/s. Maximum velocities of such values are not unlikely for energetic landslides.⁴³ However, as the validity of large maximum velocities is under debate (see also discussion at Ref. 36), the value of the drag is not further decreased. An increase in ξ results in smaller flow velocities. The deposits of scenario 3 are computed with ξ = 0.055, for which the maximum velocity of the landslide becomes approximately 50 m/s.

The value of the maximum velocity falls in agreement with values of similar events.³⁶ The simulated deposits of the scenario reach the run-out length of the deposits observed on the field [Fig. 4(c)]. Albeit, the computed deposits extend significantly beyond the lateral limits of the observed lobes in the NE-SW direction. This is possibly linked to the low yield strength of the sliding material (T₀ = 550 Pa). Indeed in scenario 4, we increase the yield strength to T₀ = 750 Pa which results in more cohesive/constrained flow deposits [Fig. 4(d)]. In this case, however, ξ has to be decreased to get an adequate run-out length (ξ = 0.045, Table I). The decrease in the drag term allows for an increase in the value of T₀. In this scenario, the simulated deposits match better the observations [Fig. 4(d)].

In scenario 4, the dynamic viscosity of the flow is also tuned to 300 Pa s, similar to the dynamic viscosity of flows with large extent as discussed in Ref. 20 (Table I). This change in the values of viscosity does not have any significant effect on the run-out length of the flow (scenario 3 compared to 4 [Figs. 4(c) and 4(d)]. This observation falls in agreement with the results of De Blasio et al., where changes in the dynamic viscosity up to two orders of magnitude in the simulations do not have significant effect on the results.²⁰ In Ref. 36, it was also observed that the best values of μ exhibit a uniform distribution.

Landslide scenario 5 simulates the landslide as a flow that exhibits behavior similar to that of a Bingham fluid but with the incorporation of a drag term. In that case, the yield strength of the landslide can be increased up to 1500 Pa for a good representation of the observations. The values of the other input parameters are kept constant (scenario 5, Table I). The computed maximum velocity of the landslide exhibits a similar trend to that of scenario 4. The simulated run-out length of the deposits falls in good agreement with the field deposits [Fig. 4(e)]. The lateral extent of the simulated deposits of scenario 5 is smaller than the extent of scenarios 3 and 4 (slightly less than scenario 4). The flow appears to be more constrained and thus in better agreement with the field observations.

From the scenarios examined, scenarios 4 and 5 form the best-fit cases for the Bingham and the Bingham with plastic terms regime (Table I). The “extended” Bingham model is possibly the best approach to simulate landslide C with VolcFlow. From a numerical point of view, increasing the basal friction of the landslide does not have any beneficial effect in the results. As the yield strength has to be reduced for the landslide to reach the desired run-out length, the sliding material spreads laterally and gives wider deposits than observed in the field (as exhibited in scenarios 3-5). However, the side scarps observed on the lower slope region constitute strong indicators of erosion and therefore dragging during the flow motion. This phenomenon renders the presence of basal friction non-negligible (with values that in reality could have been higher). For this reason, we choose the values of scenario 4 to simulate flow motion from the SUS region.

The simulated deposits of landslide A reach a run-out length similar to that of landslide C [Fig. 4(f)]. These results come in general agreement with the estimations of the extent of landslide A in the field.³⁵ However, the simulated lateral extent and run-out length of landslide A (black dashed line) are somewhat smaller than the estimated extent of the field deposits (yellow) [Fig. 4(f)]. It could be possible that the landslide in reality was more dilute than landslide C and the deposits have spread more. In the seismic profiles of the landslide deposits, the thickness of landslide A is found much smaller (∼30 m) than the thickness of landslide C.³⁵ Simulations of landslide A with a lower yield strength (T₀ = 550 Pa) and similar rheological properties, as in scenario 3, would result in deposits extending slightly more laterally.

In scenario 7, the smaller volume of the landslide (125 km³) acts as a constraining factor for the run-out length [Fig. 4(g)]. In more detail, the toe of the simulated flow comes in rest before it reaches the run-out length of landslide B (purple line) [Fig. 4(g)]. In this case, the simulated deposits spread more laterally than the field deposits. It is noted, however, that due to the limited number of available seismic profiles for the area, the estimation of the landslide deposits for landslides A and B forms an approximate interpretation of the data and in reality the deposits could have had a different extent or shape. Additional research in the area would be necessary in order to obtain more insight and advance the numerical modelling of those events.

The contours of the landslide thickness after deposition are also presented for scenarios 4-7 (Fig. 5). The thickness of the deposits for scenarios 4 and 5 appears to be fairly similar with a higher concentration of sediments closer to the landslide toe. Scenario 4 results in deposits of slightly higher thickness compared to scenario 5 [Figs. 5(a) and 5(b)]. In both cases, the maximum thickness appears to be slightly larger than 60 m which comes in good agreement with the thickness of the deposits from the seismic profiles (estimated at ca. 60 m³⁵). The thickness of the simulated deposits for landslides A and B ranges between 15 and 45 m and only rarely exceeds 45 m [Figs. 5(c) and 5(d)]. These observations also seem to match with the estimations of the real landslide thickness from the seismic profiles, which is for both cases estimated at approximately 30 m.³⁵ 

Overall, it is noted that for good agreement of the simulated run-out lengths with the field observations the values of the input parameters have to be significantly low. The low or negligible basal friction that is required for the flow to reach this run-out length may be attributed to complex processes that the code cannot capture. The same applies for the yield strength of the material. Such are the changes in the rheological composition of the mixture during motion, (especially when water is mixed with the flow) which are not uncommon in submarine sliding. Bradshaw, Tappin, and Rugg,²³ for example, have found values of T₀ = 30, 38, and 45 kPa on slope sediments, whereas they measured T₀ = 1, 10, and 20 kPa, respectively, on the deposits of the debris flow.²³ 

Other physical processes are the incorporation of water in the flow, which is only partially considered in this study. Processes that involve more uncertainty are the addition of material through erosion (which could result in volume changes during the motion) and the possible existence of an element of regression that is not visible in the current morphology. This could have an effect on the slope failure mode and possibly tsunami generation. It is also not known yet whether these events occurred independently or if they constitute composites of smaller events (i.e., comprising of smaller individual slides).

The significantly low yield strength of the modeled flow in Ref. 24 was linked with additional flow mechanisms such as remoulding and hydroplaning that can enhance the flow mobility. Hydroplaning is strongly linked with the cohesiveness and the sediment composition of the landslide.²⁰ Less cohesive flows with a sand-rich sediment composition can reach long run-outs without hydroplaning, whereas clay-rich flows tend to be more cohesive, and thus they are more likely to hydroplane.¹¹ Hence, high yield strengths may indirectly lead to greater run-out distances. However, the possibility of hydroplaning taking effect in the RBSC is not known with certainty. The erosive nature of the collapse may render hydroplaning an unlikely process for the region.

The nature of the failure mechanism and the large volume of the failed sediments raise questions for the tsunamigenic potential of the failures, as they form factors that can promote tsunami generation.⁴⁴ As shown in Ref. 36, submarine failure in the RBSC has a tsunamigenic potential that could affect the northwest Irish shoreline. As the age of most recent failure episode (landslide C) falls between 23 and 19 ka BP,³⁵ any landslide tsunamis generated in this periods would probably not have affected Ireland since the western shelf and coast were still covered by ice.⁴⁵ In more detail, the reconstruction of the different retreat stages of the BIIS and the North Sea ice cover shows that the entire country was covered by the BIIS 23 ka ago.⁴⁵ At 19 ka, some parts of the northwest coastline and shelf might have been revealed; nonetheless a significant part of the west coast was still under ice.⁴⁵ There is a high likelihood that any propagating tsunamis generated in this period would have reached an ice-covered shelf. This current hypothesis excludes the possibility of locating any tsunami deposits on the northwest Irish shoreline and therefore the validation of current and future modelling results.

It is noted that the uncertainty surrounding landslide C in the RBSC is large. Even if the above hypothesis is true, assessing the risk of submarine failure in the RBSC taking as an example the current topography of the region can be of great importance to understand the implications that a similar event could have if it occurred now. A comprehensive tsunami hazard assessment for the region is currently left for future research. In this study, we show the results of some simulations with the Nonlinear Shallow Water Equation (NSWE) solver VOLNA.⁴⁶ VOLNA gives a simplified overview of tsunami generation and propagation in the region without accounting for the complex characteristics of landslide tsunami propagation. The numerical approach of the computations (e.g., computational grid, coupling with VolcFlow) is discussed in Ref. 36. The GEBCO_08 GRID terrain model was used for the land elevation.

We present here three offshore wave gauges that measure the free surface elevation over time at water depths of 2632 m (G01), 2838 m (G02), and 1383 m (G03) [Fig. 6(a)]. The gauges record the tsunami amplitudes close to the source region and during tsunami propagation over the Rockall Trough; the tsunami simulations for scenarios 3-7 are presented [Figs. 6(a)–6(c)]. The initial water surface elevation is also given in Table II. The results show that the largest wave amplitudes are produced by landslide C. Although the initial free surface elevation between landslide A and B does not differ significantly (Table II), the discrepancy between the two scenarios increases as the tsunamis propagate [Figs. 6(a)–6(c), Table II].

Because of its lower kinetic energy, scenario 3 results in smaller tsunami amplitudes than scenarios 4 and 5 [Table II, Figs. 6(a)–6(c)]. In fact, scenarios 4 and 5 result in tsunami amplification of similar magnitude, which indicates that the rheological regime in the examined case does not affect significantly tsunami amplification. This comes in agreement with the results of rheological parametrisation studies.^21,38 The free surface elevation generated by scenarios 4 and 5 gives similar records on all the gauges with scenario 4 generating tsunamis of slightly larger amplitude than scenario 5 (Table II). A more cohesive landslide could thus result in slightly larger tsunami amplitudes. However, general conclusions about the relationship between the parameters cannot be drawn by this small sample. Furthermore, this difference is less than 5% of the tsunami amplitude of the largest scenario. Given the amount of the uncertainty that underpins the event, these discrepancies are insignificant.

The first waves that arrive in gauge 01 are the waves from scenarios 3, 4, and 5, probably due to the location of the gauge closer to the NUS-LS region [Fig. 6(d)]. The gauge starts recording the free surface elevation approximately 200 s after generation. The maximum tsunami amplitudes are 35, 34, and 32 m for scenarios 3, 4, and 5. Scenarios 6 and 7 arrive at the gauge 1 with ca. 200 s delay; the maximum tsunami amplitudes recorded are 12 and 8 m, respectively [Fig. 6(a)]. After tsunami propagation toward the Irish coastline (in a southeast direction), the tsunami amplitudes decrease {gauges 2-3, [Figs. 6(b) and 6(c)]}.

The snapshots of the free surface elevation are also presented for scenario 4 at t = 45 and t = 65 min [Figs. 7(a) and 7(b)]. The first tsunami wave crest and trough are generated between t = 50 s and t = 5 min. The tsunami amplitude is initially small (t = 50 s) but increases in amplitude and wavelength with time. As the landslide motion has a south-east direction, the first wave peak is propagating toward the Irish coastline. At the same time, the first tsunami trough also starts propagating toward an opposite northwest direction. The tsunami waves moving toward Ireland reach the eastern margin of the Rockall Trough approximately 20 min after generation.

The wave propagation is clearly confined by the morphology of the seabed over which the waves propagate. This phenomenon becomes particularly evident at t = 45 min, when the landslide tsunamis have propagated over the Irish continental shelf where the water depth decreases sharply from approximately 2000 m, at the foot of the slope, to ca. 500 m, at the shelf edge [Fig. 7(a)]. The overall form of the propagating waves exhibits a clear resemblance to the geomorphological shape of the Rockall Trough [Fig. 7(a)]. Reflection of the tsunami waves at the edge of the trough is also observed [Fig. 7(a)]. The leading tsunami waves reflect on the much shallower continental shelf, and secondary waves of smaller amplitude are formed, propagating in an opposite direction, toward the generation region [Fig. 7(a)].

The tsunami waves reach the Irish shoreline in less than an hour of propagation; the first areas to be inundated are the Mullet Peninsula, the Inishkea Islands and Achill Island, in County Mayo [Fig. 7(b)]. At t = 60 min, the leading wave crest has already inundated the northern part of the Mullet Peninsula and the first wave trough has also reached the location. Figure 7(b) shows the tsunami inundation with a smaller colour scale to increase the intensity of the inundation; it can thus be observed that the waves overtop the middle part of Mullet Peninsula. At t = 65 min, the tsunamis enter the bay in two ways: the one is by overtopping of the peninsula and the other one is by the opening between the Mullet Peninsula and Achill Island [Fig. 7(b)]. Merging of the two tsunami wave peaks and tsunami oscillations are likely to occur inside the bay at a later stage. The Inishkea Islands seem to act as barriers for the propagating waves, reducing the tsunami amplification behind them [Fig. 7(b)]. Reflection of the waves is also observed on the north part of the peninsula, the Inishkea Islands and Achill Island. After two hours of propagation, the waves have reached the Aran Islands in the south.

Overall, it is exhibited that the varying geomorphological and rheological landslide input parameters may affect significantly the depositional process of the landslide and the tsunami wave generation. Careful consideration has to be given when addressing the issue of the landslide pre-failure and post-failure conditions. Different rheological models might result in diverse landslide deposits but similar tsunami wave characteristics. On the other hand, the same rheological parameters may result in diverse wave characteristics when the geomorphological characteristics of the landslide differ substantially.

To counterbalance the lack of the water drag in the one-fluid version of VolcFLow, a drag term was introduced. In these scenarios, the peak velocities of the flow were in good agreement with the current knowledge for similar events; however, the resistance of the flow to motion was enhanced. As a result, other parameters, like the basal friction and the yield strength of the sliding material, have to be significantly reduced in order to compute the observed run-out length of the flow. Lower values of these parameters contribute to the lateral spreading of the deposits. Additional mechanisms that may enhance the flow mobility and reduction of the yield strength during the motion were not modeled by the code; these phenomena are also difficult to record in nature due to the scarce direct observations of submarine sliding.

The best-fit solutions for landslide C were obtained with the use of the Bingham rheology including a velocity dependent term (scenario 5, Table I), or the use of a rheology that incorporates the Bingham and frictional flow properties in the motion, including again a drag term (scenario 4, Table I). Comparisons of both regimes show that the Bingham rheology results in better constrained landslide deposits; however, incorporation of basal friction may be closer to the true nature of the failure. The tsunami simulations of the last episode of multiphase collapse in the RBSC with both regimes yield very similar tsunami amplitudes, albeit slightly higher in the case of a flow without friction [Figs. 6(a)–6(c)]. Landslide simulations with a frictional rheological model were also examined but were not considered advantageous for this case.

Following the best-case scenario numerical modelling of the older landslide episodes (A and B) in the south upper slope region was performed. The modeled deposits do not match as well with the observations. This can be attributed either to submarine sliding with different rheological properties (e.g., a more dilute flow) or to the current lack of knowledge regarding those events due to data limitations. Additional field and numerical research will be beneficial for a more comprehensive analysis of these episodes. In addition, a more thorough statistical approach for calibration could be used to better tune the input parameters. However, the dimension of the outputs is large since ideally both the extent curve and the thickness surface have to be compared to observations. A step in this direction has been taken recently to reduce dimension without sacrificing accuracy in the calibration.⁴⁷ The ultimate challenge would be to model, and thus calibrate, the time-varying change in the rheological composition of the mixture (and in volume assuming deposits) during motion, but these are currently intractable questions.

The landslide scenarios examined show that slope failure in the region has given rise to tsunamis of varying magnitudes. After tsunami generation, a leading tsunami crest is formed that propagates toward the Irish coastline followed by a tsunami trough. The seabed topography of the Rockall Trough strongly influences the tsunami wave propagation. This can be especially observed when the waves reach the east margin of the Trough. Reflection and refraction occur above the continental shelf as the waves pass from one medium to another. The first area to be inundated in the Irish coast is the Mullet Peninsula, in Co. Mayo. The time arrival of the first tsunami wave peak in the peninsula may be less than 50–55 min after landslide generation. At the middle part of the peninsula, tsunami overtopping is observed. The waves enter Blacksod Bay in two ways: overtopping of the land in the middle of the peninsula and propagation from the entrance of the bay.

The tsunami amplitudes and direction depend on the scenario examined (volume, maximum thickness of the landslide, region of collapse). Focusing on the three major episodes of multiphase collapse attributed to the RBSC, the most voluminous event (landslide C) has the largest tsunamigenic potential. Gauges close to the generation region measure a maximum tsunami amplification of approximately 40 m that decays significantly during propagation. Landslide-induced tsunamis from the SUS region are smaller in amplitude. The second episode of collapse (landslide B) holds the smaller tsunamigenic potential of the three.

The results of this study constitute an attempt to study numerically the different episodes of collapse in the RBSC under a rheological framework and offer a view on subsequent tsunami generation. Future research could be devoted to a more sophisticated modelling of the phenomena. This could be approached by employing a fully coupled numerical code. In addition, to adequately assess the tsunami propagation and predicted run-up heights on the coast, the use of more advanced dispersive solvers can improve the modelling in some settings. Finally, a slope stability assessment at the source and a different meshing strategy in the modelling with mesh refinement on the coastal areas (where impact would be better estimated as a result) would be of great value for a comprehensive hazard assessment in the area.

The authors are thankful to Dr. Kelfoun for kindly providing with the VolcFlow code for the landslide modelling. This research was supported by the European Union’s Seventh Framework Programme for research, technological development, and demonstration under the grant agreement ASTARTE No. 603839. D. M. Salmanidou and S. Guillas also gratefully acknowledge partial funding from the M2D network (Grant No. EP/P016774/1, Award No. M2DPP017, Project title: “Potential large tsunami hazards associated with landslide failure along the West coast of India: from uncertainties to planning decisions”).