Previously undiscovered landslide deposits in Harrison Lake, British Columbia, Canada

A bathymetric survey of Harrison Lake in southwest British Columbia revealed deposits of three large landslides on the lake floor. The blocky and flow-like surface morphology of the deposits suggests rapid emplacement from subaerial sources. The multibeam survey, together with a subbottom acoustic survey, allowed us to estimate deposit volumes of 2.4 Mm3, 1.3 Mm3, and 0.2 Mm3 for the Mount Douglas, Mount Breakenridge, and Silver Mountain landslides, respectively. The large volumes and inferred rapid emplacement of the Mount Douglas and Mount Breakenridge landslides suggest they were tsunamigenic. Because people live along the shoreline of Harrison Lake, our discovery and characterization of these landslide deposits and their tsunami-generating potential form an important foundation for further landslide-tsunami hazard analysis in the region.


Introduction
Fjord coastlines, notably those in Canada, Alaska, Greenland, Norway, Chile, and New Zealand, are vulnerable to landslidegenerated waves (Roberts et al. 2014;Gauthier et al. 2018;Higman et al. 2018). Additionally, there is a growing awareness of the hazard posed by landslides into lakes bordered by steep, potentially unstable slopes in mountains worldwide (e.g., Roberts et al. 2014;Haeberli et al. 2017).
Landslide-generated waves can be much larger than earthquake-triggered tsunamis, at least near their sources (Tappin et al. 2001;Graziani et al. 2006;Omira et al. 2019) and are common worldwide (Roberts et al. 2014). Historic events exemplify their potential magnitude. The 1934 Tafjord (Norway) rockslide generated a displacement wave that reached up to 63 m above sea level (asl), and the 1936 Lake Loen (Norway) rockslide produced a wave with a 74 m run-up (Grimstad 2005;Oppikofer et al. 2009). The 1963 Vajont rockslide (Italy) occurred in a reservoir and generated a displacement wave, which overtopped the reservoir damn by 245 m, causing an estimated 2000 casualties (Müller 1964). The 1958 Lituya Bay (Alaska) rockslide produced a wave that reached up to 517 m above sea level (Miller 1960), and the 2015 Taan Fiord (Alaska) rockslide generated a 193-m-high wave at the head of the inlet Haeussler et al. 2018;Higman et al. 2018). In 2007, a 3 Mm 3 rockslide entered Chehalis Lake, 80 km east of Vancouver, British Columbia, generating a displacement wave that reached up to 38 m above the lakeshore and ran down the 8-km length of the lake (Brideau et al. 2012;Roberts et al. 2013). The Chehalis lake event, in part, motivated our research on the exposure risk of large subaerial landslides in British Columbian lakes. Given the apparent increase in large landslides in British Columbia (BC) in recent years, due in part to loss of glacier ice and permafrost thaw (Geertsema et al. 2006;Cloutier et al. 2016;Coe et al. 2018;Hibert et al. 2019), the potential exists for more frequent tsunamigenic landslides in the future (Roberts et al. 2014;Higman et al. 2018).
In this paper, we document evidence of three previously unreported landslide deposits in Harrison Lake, southwestern BC ( Fig. 1), that we discovered during a bathymetric survey of the lake in 2017 and 2018. Using a bathymetric digital elevation model (DEM) and seismic imaging of the lake sediments, we characterize the morphology of the landslide deposits, infer emplacement mechanisms, and apply a statistical approach to estimate deposit volumes. Our findings provide a foundation for further study of hazard and risk of landslide-generated waves in Harrison Lake, and other lakes in BC. Our research is motivated by the ongoing mandated commitment of the British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural Development (FLNRORD) to conduct research on natural hazards in British Columbia.

Setting
Harrison Lake-90 km east of Vancouver, British Columbia and < 10 km east of Chehalis Lake-is the easternmost and largest (218 km 2 ) of a series of northerly trending fiords and fiord-lakes between the Salish Sea and Fraser Canyon (Fig. 1). It has a length of~60 km, a maximum width of 8 km, and a maximum depth of 280 m (obtained from our survey). The surface elevation of Harrison Lake is 10 m asl, whereas the sleep slopes surrounding the lake rise to elevations above 2000 m asl. Lillooet River flows into its north end of the lake and Harrison River, a tributary of Fraser River, flows out to the south at Harrison Hot Springs, a community of 1500 people. The lake has a low specific sediment yield (0.102 mg km −2 day −1 ; Slaymaker and Menonous 2000); thus, Holocene landslide deposits are well exposed on the lake floor adjacent to the steep valley walls. Fluvial and deltaic deposits are present where streams enter the lake, and lower slopes are mantled by colluvium and thin till.
The bedrock geology of Harrison Valley is a product of Mesozoic accretion and subsequent compression and plutonism of the southeastern Coast Mountains (Monger and Price 2000). Faults separate Jurassic and Cretaceous volcanic, marine sedimentary, and granitic rocks west of the lake from metamorphic rocks of varying grades and granites to its east (Gabites 1985;Journeay and Friedman 1993;Mahoney et al. 1995;Brown et al. 2000;Cui et al. 2017). The three landslide deposits that we describe in this paper are sourced in Slollicum Schist (Cui et al. 2017).

Bathymetric and subbottom data acquisition
We conducted a multibeam echosounder (MBES) survey of Harrison Lake from 8 to 17 July 2017 and 7 to 13 July 2018. In 2017, we used a 400 kHz Reson® 7125 MBES with a maximum 165°acrosstrack swath angle and a ping rate of 2 Hz (deep water) to 50 Hz

Original Paper
Landslides 18 & (2021) (2021) 18:529-538 (shallow water). We georeferenced and corrected for vessel motion using an Applanix PosMV V5 system consisting of an inertial measurement unit (IMU) and Trimble Global Navigation Satellite System (GNSS) receivers. In 2018, we used a Norbit WBMS MBES with a maximum 210°across-track swath angle, off-nadir capability, and ping rate of 2-50 Hz. We georeferenced and corrected for vessel motion using an integrated IMU and GNSS (Applanix Ocean Master). We applied measurements of offsets between the instrument components in the acquisition software PDS 2000, and used a sound velocity probe to collect 35 vertical profiles (casts) to correct the acoustic signal travel time of the MBES system. We used Qimera software to process and clean apply the raw data, gridded the point data to a 1-m raster surface, and interpolated data gaps via adaptive triangulation in Quick Terrain Modeler software.
We acquired acoustic reflection profiles of the subbottom deposits between 11 and 13 June 2018 using an Edgetech 3200 CHIRP acoustic subbottom profiler operated at 2-16 kHz. The Edgetech navigation system malfunctioned during the survey; thus, we had to georeference the acoustic reflection profiles based on the alignment of timestamps with the Lowrance GPS onboard the survey vessel, which resulted in a positional accuracy of 5 m We processed and interpreted the collected data using Edgetech's Discover proprietary software and an assumed sound velocity of 1500 ms −1 to convert two-way travel time to depth. No cores were collected for this study; therefore, our inferences of the characteristics of the  subsurface sediments, as well as their seismic velocities, are based solely on interpretation of the seismic imagery.

Landslide dimensions
We delimited the margins of landslide deposits and generated polygons on the basis of differences in their ridged, raised, and blocky form and adjacent smooth lake floor. We used the subbottom acoustic reflection data to corroborate the deposit margins by identifying the boundaries between the undisturbed, laminated, acoustically transparent lacustrine facies and the hummocky, acoustically impenetrable reflectors of the coarse landslide deposits. After draping landslide polygons over the bathymetric DEM, we calculated travel angles, runout lengths, and subaqueous deposit extents in a GIS.
The common approach to quantify landslide volumes involves differencing measured pre-and post-landslide surfaces (Cruden and Varnes 1996;Gauthier et al. 2018;Haeussler et al. 2018). In our case, the pre-landslide surfaces are unknown; therefore, we estimated landslide deposit volumes by modeling pre-landslide bathymetry and differenced those surfaces from the measured postlandslide bathymetric surfaces, an approach similar to that used in previous studies (Iverson et al. 2015;Sandøy et al. 2017;Cui et al. 2018).
We reconstructed the pre-landslide bathymetry in four stages. First, we determined the slope and curvature of the lake floor surrounding the deposits based on a polynomial regression analysis of longitudinal (downslope from lake shore to depth) and lateral (contour parallel to shore) profiles. Polynomial models of different orders were tested, as local surface representations of curvature require higher orders to reflect bathymetric complexity (Schmidt et al. 2003). The best-fit models had R 2 values of 0.95 to 0.99, and an F test P value lower than the 0.05 significance level. Second, we clipped and removed the data contained within the landslide polygon. Third, we generated a mesh of lateral and longitudinal profiles over the clipped area and assigned depth values based on the slope and curvature regression models. The generated depth values were constrained to the depth range in the bathymetry, i.e., no deeper than maximum depth in the surrounding area could be created. Finally, we merged the mesh with the manipulated clipped point cloud, regridded the data to 1 m using the local polynomial regression interpolation method and applied a smoothing filter in a GIS. We then performed DEM differencing using the Geomorphic Change Detection (GCD) Tool for ArcGIS, with a 0.2 m threshold (Wheaton et al. 2010). Elevation points from the measured 2017/2018 bathymetry that were above the prelandslide surface were considered deposition and points below were considered erosion. The "total volume of mobilized sediments" was based on all elevation changes, both positive and negative, and the net volume was the deposited volume minus the eroded volume, therefore the total volume that entered the lake. The GCD also provided an estimate of deposit thickness calculated from the average pixel-to-pixel differences.
Several assumptions are implicit in our analysis of landslide deposit volume. First, we assume that the lake bathymetry beneath each landslide deposit is the same as the surrounding terrain, and thus ignore the possibility of constructional (e.g., fan) or erosional (e.g., channel or gully) forms that may have been buried by the landslide deposits. Second, we assume that each deposit is a product of a single slope failure and not of multiple events. Third, we assume that sedimentation is negligible, which may not be the case in some areas, notably near Silver Mountain, which is near Silver River and Bear Creek. Lastly, we assume there was no prior erosion where there has been subsequent deposition.

Landslide detection
The multibeam data show that the landslide deposits are elevated above the adjacent lake floor and have morphological characteristics typical of subaqueous landslide deposits (Clare et al. 2019), including arcuate transverse ridges, hummocks, and angular blocks (Figs. 2, 3, 4, and 5). We interpret these deposits to be subaerially, as opposed to subaqueously, sourced due to the prevalence of angular blocky debris and upslope instabilities. This interpretation is supported by comparisons to analogous events (Hermanns et al. 2014b;Haeussler et al. 2018), as well as the absence of potential subaqueous source areas in the bathymetry.
The seismic data corroborate our interpretation of these features and allowed us to delineate landslide margins. The nonlandslide lacustrine sediments exhibit horizontally stratified layers with high acoustic transparency, whereas the landslide deposits are characterized by hyperbolic and chaotic reflectors, indicative of large blocks in the subsurface (Fig. 5). High acoustic impedance of the landslide material limited the penetration of energy into the subsurface and prevented us from establishing the thickness of the deposits based on the seismic data.

Landslide descriptions
We describe the subaqueous components of the three landslides based on interpretation of the MBES bathymetry and subbottom survey data. The source-area boundaries of the landslides are not precisely resolved and await further field and Lidar investigation and interpretation.

Mount Douglas
The Mount Douglas deposit is the northernmost and largest of the three deposits, with an area of 0.8 km 2 (80 ha) (Figs. 1 and 2). It has a runout distance of 1.5 km, a maximum width of 1.2 km, and estimated volume of 2.4 Mm 3 (Tables 1 and 2). Our analysis suggests the deposit has an average thickness of 4 m and a maximum thickness of 16 m ( Table 1). The source slope extends upward to an elevation of 1300 m asl with an average inclination of 36 o ; the underwater portion of the slope is inclined 22°( Table 1). The subaqueous travel angle is 9° (Table 1, Fig. 2D).
The deposit is a dense aggregate of large blocks with some prominent isolated blocks at its margins (Fig. 2B). Transverse ridges with an average relief of 1.5 m are located in the upper flow path and are flanked by two lobes of sediment derived from neighboring gullies ( Fig. 2A). Conspicuous distal arcuate ridges (average height = 5 m) and shear zones are present along the south margin of the deposit (Fig. 2C). The subbottom profile (Figs. 4A and 5) indicate hard reflectors transitioning to laminated lake sediments and show the hummocky nature of the deposit.

Mount Breakenridge
The Mount Breakenridge deposit is located 9 km south of Mount Douglas (Figs. 1 and 3). It has an area of 0.5 km 2 (50 ha), a runout distance of 1.3 km, a maximum width of 1.0 km, and an estimated volume of 1.3 Mm 3 . Our analysis yields an average deposit thickness of 4 m and a maximum thickness of 27 m, with hummocky hard reflectors transitioning to laminated lake sediments in the subbottom profile (Fig. 5). The source slope extends to a maximum height of 1350 m asl with an average inclination of slope 35°; the subaqueous portion of the slope is inclined 23°( Table 1). The subaqueous travel angle is 8° (Table 1, Fig. 3D). Our MBES and subbottom surveys document bulging of the submerged toe of the slope at Mount Breakenridge, 500 m south of the deposit (Fig. 6).
Isolated blocks~9 m high and 25 m wide are present on the steep subaqueous slope above the main depositional area (Fig. 3A). The thickest portion of deposit is located at the base of this steep slope (Fig. 3B) and supports lenticular ridges up to 30 m long, 5 m wide, and 3 m high. The ridges have flow-parallel orientations, but individual fragments are transverse to flow. More distally, within the main flow path, angular flat-topped blocks up to 30 m long and 4 m high have a stepped/terraced form (Fig. 3C). Beyond this zone is an outermost zone comprising smaller blocks.

Silver Mountain
The Silver Mountain deposit (Figs. 4 and 5C) is 17 km south of Mount Breakenridge. British Columbia archival photographs illustrate a fresh non-vegetated slope at this site in the late 1800s (City of Vancouver Archives item CVA:1376-375.38; Bailey and Neelands 1889). The deposit has an area of 0.17 km 2 (17.9 ha), a runout distance of 0.8 km, a maximum width of 0.9 km, and an estimated volume of 0.2 Mm 3 (Tables 1 and 2). Its average and maximum thicknesses are, respectively, 1 m and 12 m, and thus, it is much thinner than the other deposits. The source slope is lower (< 850 m . The red lines delineate the deposits margins, and the blue line the location of the subbottom data transect, the white line the centre flow path profile for inset D. Inset A depicts small transverse ridges, inset B prominent isolated blocks, inset C arcuate compressional ridges, D a central amorphous zone, E a zone of translational blocks, F a lobe of smaller blocks. Inset D illustrates the elevation profile over the central flow path asl) and substantially less steep (23°) than those at the other two sites. The subaqueous portion of the slope has an average angle of 26°. The subaqueous travel angle is 13° (Table 1, Fig. 4D).
The morphology of the Silver Mountain deposit is much different from that of the other deposits. The deposit is wide and has a limited runout, apart from one distinct lobe to the south. The deposit consists primarily of low-relief transverse hummocks, which indicates that it lacks the blocky material observed in the other deposits (Fig. 4A). Small isolated angular blocks are present only in the protruding lobe at the south margin of the deposit (Fig.  4C); the remainder of the deposit has an amorphous form with muted microtopography (Fig. 4B). The subottom profile (Fig. 5) also indicates a thinner deposit with an acoustically transparent character, suggesting a paucity of rocky rubble.

Discussion
Landslide deposit morphologies and emplacement mechanisms Documentation of deposit morphologies and comparisons to previously studied comparable events allow us to infer landslide transport and emplacement mechanisms. The deposit morphologies for all three events reflect rapid runout and deposition from a subaerial source.
The Mount Douglas deposit is morphologically similar to deposits of other large subaerial landslides that have entered fjords (e.g., Hermanns et al. 2014b;Haeussler et al. 2018). The long runout, ridge morphology, large blocks, and proximity to steep entered Harrison Lake and likely eroded and incorporated additional material sourced from gullies directly north and south of the deposit. Deposition in this area was limited due to the high velocity of the mass on its impact with the lake. The debris sheet spread as it crossed the lower slope of the lake floor. The velocity of the flowing debris mass was highest along its centerline, and much of the debris moved to lateral and distal areas (e.g., Shea and van Wyk de Vries 2008), leaving a muted amorphous central zone ( Fig. 2A; Geertsema et al. 2018). The arcuate distal ridges formed due to compression as the debris mass slowed (de Blasio and Mazzanti 2010). Similar features are common in both historical subaerial and subaqueous rockslide and rock avalanche deposits, for example the deposit of the 2007 rockslide in Aysén Fiord (Hermanns et al. 2014b) and the 1881 subaerial rock avalanche at Elm, Switzerland (de Blasio and Mazzanti 2010). Another morphological indicator of rapid flow behavior is the transverse ridges in the upper flow path (Baloga and Bruno 2005).
The presence of distal out-runner blocks is evidence of translational sliding within the flowing debris mass (Fig. 2, A3). Outrunner blocks are large fragments of source material that travel beyond the main debris mass due to their high momentum and are found outside the main debris margins of many subaerial and subaqueous rockslide and rock avalanche deposits (de Blasio and Mazzanti 2010). The blocks in Fig. 2 (A3) are similar to those at the distal margin of the 2015 Taan Fiord landslide .
The Mount Breakenridge deposit morphology differs from that at Mount Douglas, suggesting a different emplacement mechanism. The main area of accumulation is proximal to the source slope (Fig. 2, B2). The shorter runout length and the position of the main accumulation proximal to the source, rather than at the distal margins, indicate lower mobility. The transport velocity was insufficient to sustain movement, and thus the debris slowed and became limited to sliding motion (Hungr et al. 2014).

Original Paper
The surface morphology of the Mount Breakenridge deposit is ridged, with flow-parallel orientations and transverse fragments within ridges. Slope changes and irregular micro-topographic variability of the lake floor likely affected the momentum of the flowing debris, creating these ridges (Baloga and Bruno 2005). Longitudinal ridges are found in landslide deposits with nonuniform grain size distributions and have been attributed to flow changes over a high-friction substrate (Dufresne and Davies 2009). The smaller transverse ridge fragments may indicate a change in velocity attributable flow instability, clast-size differences, surface irregularities, or pulses in the release of source material (Baloga and Bruno 2005).
Translational sliding is evidenced by the presence of the flat-topped blocks with a stepped form in the distal part of the Mount Breakenridge deposit. The blocks are not of lacustrine origin as they have high acoustic impedance (Fig. 5) and thus are probably rock. We suggest that large slabs of schistose bedrock slid on lakebed sediments and became imbricated in the distal part of the runout zone. The longer runout of this nearly intact source material is attributed to its high momentum (de Blasio and Mazzanti 2010). The Mount Douglas and Mount Breakenridge deposits are derived from the same geological unit, but have different morphologies, which implies different emplacement mechanisms. The Mount Douglas deposit has a morphology indicative of flowing behavior that we ascribe to incorporation of subaqueous slope material and subsequent bulking, whereas the Mount Breakenridge deposit was emplaced mainly by sliding. In the Mount Douglas case, the greater amount of granular material, derived either from the source or by entrainment from the lake floor, would have led to motion closer to that of a rock avalanche, in which flow behavior is dominant, than a rockslide (Shea and van Wyk de Vries 2008). At Mount Breakenridge, the prevalence of intact source material and limited run-out indicate sliding behavior.
The Silver Mountain deposit is the smallest of three landslide deposits documented in this study and lacks the morphological features characteristic of the other deposits. The archival photograph of the source slope (Bailey and Neelands 1889) suggests that it was emplaced by a debris avalanche involving till, colluvium, and perhaps weathered bedrock. It shows a shallow scar that widens downslope. The thinness of the deposit, the hummocky morphology, the large width-to-runout ratio, and small volume support a debris avalanche interpretation (Hungr et al. 2014).

Implications for landslide-generated waves
Comparison of the volume and morphology of the Harrison Lake landslide deposits with those of well-documented historic events such as the nearby 2007 Chehalis Lake rock avalanche (Brideau et al. 2012;Roberts et al. 2013) suggests that the Mount Douglas and Mount Breakenridge landslides may have produced sizable displacement waves within the lake. Hermanns et al. (2014a) cataloged landslide-generated waves in Norway and found that landslides with volumes greater than about 0.5 Mm 3 could generate displacement waves with run-ups greater than 20 m and measurable effects more than 10 km from the source. The Mount Douglas and Mount Breakenridge landslide deposits exceed this threshold volume.
Numerous landslide-generated waves have been documented around the world (Roberts et al. 2014 (Grimstad 2005;Oppikofer et al. 2009;Brideau et al. 2012;Roberts et al. 2013). The tsunami potential of the much smaller Silver Mountain debris avalanche (~0.25 Mm 3 ) is inconclusive. However, even small, fast-moving events can produce localized sizable waves with run-ups that do not always scale simply with failure volume. For instance, a 0.15 Mm 3 volcanic debris slide (Wagner et al. 2003) and a 0.1 Mm 3 rock avalanche (Plafker and Eyzaguirre 1979) produced proximal waves up to, respectively, 15 m and 30 m high.
The occurrence of two, possibly three, wave-generating landslides on the east flank of Harrison Lake suggests that the stability of the source slopes should be better analyzed and characterized. Field observations of flexural topples and slope parallel fractures on Mount Breakenridge (Nichol et al. 2002) and our documentation of a bulging toe slope beneath the lake surface (Fig. 6) suggest a long history of slope deformation. Table 1 Metrics of landslide deposits measured from MBES bathymetry, with the exception of the above water hillslope angle measured from Google Earth. Thickness measurements obtained from Geomorphic Change Detection using 1.5 m raster pre-landslide surfaces and 1. Average thickness (from GCD) (m) 4 ± 5% 4 ± 6% 1 ± 16% Maximum thickness (from GCD) (m) 16 27 12 Volume of mobilized sediments (total volume) (m 3 ) 2,700,000 ± 150,000 1,500,000 m 3 ± 94,000 m 3 310,000 ± 58,000 Landslide volume (net volume) (m 3 ) 2,400,000 ± 140,000 1,300,000 ± 80,000 220,000 ± 81,000 Original Paper

Conclusion
We conducted multibeam and seismic surveys to map the bathymetry and sediment characteristics of Harrison Lake in 2017 and 2018. During these surveys, we found three large, previously undocumented landslide deposits on the east side of the lake at Mount Douglas, Mount Breakenridge, and Silver Mountain. The Mount Douglas and Mount Breakenridge deposits are the result of, respectively, a large subaerial rock avalanche and a rockslide that entered the lake at high velocities and ran out more than 1 km across the lake floor. Based on other studies, we conclude that these two events were likely generated large displacement waves. The much smaller Silver Mountain deposit was a product of a debris avalanche and likely produced a small displacement wave. Our bathymetric data also reveal slope bulging at the base of Mount Breakenridge, indicating that the slope continues to slowly deform and could be the site of a future landslide. We further recommend dynamic modeling of the past landslides together with wave modeling based on our calculated landslide volumes and surveyed lake bathymetry.

Funding information
Funding for student salaries, equipment rentals, and field expenses was provided from a FLNRORD IO grant to Geertsema. The 2017 loan of the MBES was made possible by a NSERC grant to Jeremy Venditti of Simon Fraser University through the River Dynamics Lab.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.