We describe and model the evolution of a recent landslide and outburst flood in the southern Coast Mountains, British Columbia, Canada. About 18 Mm³ of rock descended 1000 m from a steep valley wall and traveled across the toe of a glacier before entering a 0.6 km² glacier lake and producing a >100-m high wave. Water overtopped the lake outlet and scoured a 10-km long channel before depositing debris on a 2 km² fan below the lake outlet. Floodwater, organic detritus, and fine sediment entered a fjord where it produced a 70-km long turbidity current and altered turbidity, water temperature, and water chemistry for weeks. The outburst flood destroyed forest and culturally significant salmon spawning and rearing habitat. Physically based models of the landslide, the displacement wave, and the flood provide real-time simulations of the event and can improve understanding of similar hazard cascades and the risk they pose.

Airborne laser altimetry

We performed three airborne laser altimeter surveys over the study area using a Riegl Q780 full waveform scanner (1064 nm). As part of a larger forest inventory survey, a commercial party acquired lidar data on 28 August 2014 for the lower reaches of Elliot Creek (up to about 2 km upvalley from the mouth of the creek). These data were combined with pre-event lidar data obtained during a survey of the Homathko Icefield on 11 October 2019. Finally, post-event lidar data were acquired on 23 December 2020 with the Riegl Q780 scanner mounted on a helicopter. Average point density for the three surveys are 15, 2.2, and 25 laser shots points/m⁻², respectively. We classified the lidar point data into ground and non-ground laser returns and gridded the latter into 1 m bare earth Geotiffs. We co-registered the two pre-event lidar grids and subsequently co-registered the pre-event data against the post-event lidar product. Due in part to dense vegetation cover and the low point density of the October 2019 lidar survey, we ascribe an uncertainty within vegetated terrain of ± 0.69 m (1σ). In contrast, vertical and horizontal uncertainty over sparsely vegetated and bare surfaces is closer to ± 0.25 m. Although our true positional errors are likely lower given the large sample size, our surveyed area is too small to correct our sample size for spatial autocorrelation. We thus conservatively use ± 0.69 m to quantify the vertical and horizontal uncertainty of these surveys.

Landslide and tsunami modeling

We used a hydrodynamic model (FLOW-3D HYDRO) to simulate the whole process of landslide generation and resulting tsunami propagation in Elliot Lake.

FLOW-3D HYDRO is a CFD software that solves the hydrostatic Reynolds-averaged Navier-Stokes equations to obtain transient, three-dimensional solutions to multi-scale, multi-physics flow problems. A key feature of this numerical model is the implementation of the Volume of Fluid (VOF) method for simulating free surfaces . The VOF method is a numerical technique used to track the locations and movement of complex free surfaces and to apply proper dynamic boundary conditions to those surfaces. The model uses a structured computational mesh that is composed of rectangular elements defined by a set of planes perpendicular to each of the coordinate axes.

Critical inputs to FLOW-3D are topographic data and the modelling approach to simulate the rock slide. We imported the topographic lidar data from 11 October 2019 (pre-slide event) into the program as a raster (Figure S11). Given the modelling area and computer capacity, we used a 10 m x 10 m grid size for the model. The water surface elevation of Elliot Lake was set to 727.5 m asl, which corresponds to an approximate lake surface area of 59 ha. We used a 2D radial basis function interpolant to estimate the lake bathymetry using the adjacent side moraines and slopes to guide the interpolation.

The rock slide was simulated as a clear-water flow with a total volume of 13.3 Mm³, which is equal to the estimate of the deposit in the lake. The model outflow boundary is approximately 300 m south of Elliot Lake and was set to an outflow-type boundary in FLOW-3D. We used the renormalized group (RNG) model to describe the turbulent flow and set a surface roughness height to 0.01 m for the entire domain.

The time series data for the resulting tsunami as it propagates away from the impact zone indicates there are multiple reflected waves that continue to oscillate within the lake after the leading wave has overtopped the dam. These reflected waves have a much smaller maximum amplitude and a significantly longer wave period, which translates to a slower moving wave.  The time series data also indicates that the maximum amplitude of the leading wave decreases with time until P5 (x = 1965 m). 

Outburst flood modeling

There were no direct measurements or observations of the outburst flood, and so our estimates of its properties came from a reconstruction approach, sometimes referred to as palaeoflood analysis. We employed Delft3D model software (version 4.02.03) in this study because it can accommodate a wetting front (i.e. propagation of a wave over an initially dry bed) and hydraulic jumps (rapid transitions between sub- and super-critical flow regimes). It also can simultaneously consider both hydraulic variables and interactions between the flow and the bed, enabling both sediment transport and iterative bed elevation updating (i.e. erosion and deposition at each model time step). Specifically, the physically based, fully non-linear, Navier-Stokes equations were solved based on shallow water assumptions and computed on nodes of a curvilinear computational mesh that was created based on high-resolution (2 m), pre-flood topographic data acquired during an airborne lidar survey. Design of the model domain (i.e. the extent and coverage of this curvilinear mesh) was guided by the mapped wash limits.

In this study, numerical simulations were used to calculate hydraulic parameters including depth-averaged velocity, flow depth, Froude number (Fr), and bed shear stress (ꞇ) for every node of the computational mesh and for every model time step, which in order to conserve mass and momentum, was necessarily small (0.0001 mins).

The outburst flood model was parameterised assuming a spatially uniform or ‘global’ roughness of Manning’s n of 0.05, which in lieu of spatially distributed grain-size information, we consider to be representative of glaciated mountain valley gravel-bed rivers. Similarly, in lieu of any knowledge of distributed sediment depth and grain size distributions, we assumed a spatially uniform or ‘global’ sediment depth of 20 m, a specific density of 2650 k gm^-3, and a D50 grain size of 2000 µm, which we consider to be representative of the alpine glacial valley sediment in this region of western Canada. Although the model incorporates the effect of sediment in the water column on water density, it is a Newtonian fluid model only, so there is no consideration of sediment concentration effects on flow behaviour. Additionally, only vertical exchanges of sediment with the bed are considered; there is no formulation or representation of mass failures such as bank collapses.

We used two conditions to introduce water to our outburst flood model. The first is the height of the mapped wash limits above the valley floor at the moraine dam, which indicate peak flood stage at our outburst flood model upstream boundary. The second is the output hydrograph of the tsunami model, which indicates the rate of rise to peak discharge, as well as the rate of the falling limb, at the moraine dam. Our outburst flood model was therefore forced by a time-transgressive water level, reaching peak stage just 30 seconds after the entry of the landslide into the lake and then exponentially falling to valley floor level after 3 minutes. The downstream boundary of our outburst flood model was sea level and was parameterised as such; i.e. water level equal to 0 m a.s.l. for the duration of the flood. A ‘spin-up’ model was used as input to the full model runs. The spin up model was a 48-hour run using a base flow of 10 m³ s^-1 , which we regard as a realistic background discharge for Elliott Creek at the time of the landslide, given scaling of comparable discharges in the region based on catchment size. This model served to ‘pre-condition’ or re-distribute the mobile sediment realistically across the model domain.

Fjord turbidity and water temperature analysis

Temperature, salinity, pressure (CTD), and oxygen data have been collected at eight evenly spaced stations that span the length of Bute Inlet since 1951 by the University of British Columbia, Fisheries and Oceans Canada, and the Hakai Institute. Following Jackson et al. (2021), we removed the seasonal cycle from temperature, salinity, and oxygen data by using monthly averaged data from 1951 to 2010. Once the seasonal cycle was removed, the time series could be examined in the context of stochastic, annual, and decadal change. For the analysis presented in Figure 4, data were collected 37 days before (22 October 2020), four days after (2 December 2020), and 16 days after (14 December 2020) the Elliott Creek landslide. 

Turbidity data and calibration  

Beginning in 2018, Seapoint turbidity meters were attached to the CTD systems to investigate particle dynamics within Bute Inlet. These sensors have a light source in the near-infrared (880 nm), were factory calibrated with Formazin (output in FTU), and were rated to have a linear response (+/- 2%) to particle concentrations spanning 0 to 1250 FTU. Unfortunately, a different sensor was utilized on the two surveys following the landslide, and this sensor showed an approximately 2 FTU greater offset than the sensor utilized prior to the landslide. A sensor-specific offset, rather than increases in particle concentrations between surveys, was proven by comparing data from the two sensors in Bute Inlet with data collected at stations outside Bute Inlet, both of which showed a comparable offset over a wide range of conditions. As a result, corrections were performed for this offset to make the data from both sensors comparable. First, we selected a station that consistently showed the lowest profile turbidity magnitudes over the time series as a reference station (BU2, 658 m deep and 54 km from the head of the fjord). Second, for each survey, we averaged the 10 minimum turbidity values from the profile collected at this reference station to represent a “particle-free” instrument blank value. Data were averaged over the 10 minimum values to limit the influence of instrument noise, and the standard deviations of these averages were low (0.002–0.006 FTU). Third, this derived instrument blank value was then subtracted from data for each station collected during the survey. At times, this method resulted in slight negative values; these were set to zero. Care should be taken when interpreting these data, as it is unlikely that true “particle-free” blank offset values could be determined in dynamic coastal waters with high and variable particle concentrations. However, the calculated offset values were consistent over time in comparison to stations outside Bute Inlet, suggesting a lack of environmental influence and adding confidence to the method.

The 6 animations of various aspects of the landslide, flood, and tsunami (impulse wave) are in .avi or .mp4 format.  These are meant to be visualizations to help understand the dynamics of the cascading hazards.

The lidar difference data is in .asc format, readable by open source GIS software.

Fjord turbidity and CTD (temperature, salinity, and pressure) data are available in .csv and.xlsx formats.