Numerical modeling of a potential landslide-generated tsunami in the southern Strait of Georgia

We report results of numerical simulations of a potential subaerial landslide on the coast of Orcas Island and the resultant tsunami waves in the southern Strait of Georgia near the US/Canada border. A likely trigger is strong ground shaking during large earthquakes on the nearby Holocene active Skipjack Island fault zone. For a worst-case scenario, we assume a 0.17 km3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{km}}^3$$\end{document} rigid subaerial failure on the steep northeast coast of the island, spanning the ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document} 5 km between previous landslide deposits on the adjacent seafloor. The landslide motion and resulting tsunami generation are modeled using the three-dimensional (3D) non-hydrostatics physics-based NHWAVE model. The simulated failure moves downslope with a peak velocity of 13.64ṁ/s and travels 732 m before coming to rest after 85 s in 75-m water depth. Tsunami propagation is then continued using the 2D fully nonlinear and dispersive Boussinesq wave model FUNWAVE-TVD in a succession of layered and nested grids. The modeling reveals susceptible locations, particularly as waves will arrive with little or no warning. In the near-source region, modeled waves have peak amplitudes of 15–20 m, current speeds of up to 10 m/s, and runup of up to 30 m. Smaller, but significant, wave amplitudes and runup occur throughout the region surrounding Orcas Island. In the tsunami propagation direction, runup reaches 7.5 m at Neptune Beach near Lummi Bay. Both initial and reflected waves cause significant runup (> 1.5 m) along much of the shoreline between Point Roberts and Lummi Bay.


Introduction
The coastline of British Columbia (B.C.) is at risk from tsunami waves generated by a range of different sources, with most studies focusing on seismically generated tsunamis that have the potential to cause the most widespread damage (e.g., Cherniawsky et al.

3
2007; Fine et al. 2008;Gao et al. 2018). However, it is also critical to assess and mitigate the tsunami risk posed by subaerial and underwater landslides.
Several events over the past century have demonstrated that mass movements such as submarine slumps and subaerial rockfalls play a far more significant role in tsunami generation than previously thought. Many "tsunami earthquakes" are followed by landslides, which locally increase the tsunami waves generated by tectonic displacements (Iwasaki et al. 1996). Examples include the 1922 Flores Island Gica 1996), 1998 Papua New Guinea (Imamura et al. 2001(Imamura et al. ), 1999 Kocaeli, Turkey (Altinok et al. 1999), and 2011 Tohoku, Japan (Grilli et al. 2013) events. Landslide-generated tsunamis (LGTs) are far more confined than those induced by earthquakes, and they can cause extreme coastal runup and associated damage, particularly when the wave energy is retained by inlets or semi-enclosed embayments (Rabinovich et al. 2003). Tsunami runup resulting from subaerial landslides can be particularly high, as evidenced by several historical events with runup exceeding 100 m (e.g., Higman et al. 2018).
Landslides and triggered waves tend to recur regularly in some areas. After the 525 m runup tsunami in Lituya Bay, Alaska in 1958, field surveys indicated that LGTs had previously occurred in the bay in 1853-54 (120 m), 1874 (24 m), 1899 (60 m), and 1936 (150 m) (Miller 1960;Rabinovich et al. 2003). Other recurring LGT sites along the west coast of North America include Kitimat Arm, B.C. (Johns et al. 1985), Yakutat, Russell Fjord, Skagway Harbor (all in Southeast Alaska) (Kulikov et al. 1996;Lander and National Geophysical Data Center 1996), and Puget Sound (Washington State) (Baum et al. 2005).
Local subaerial (e.g., Knight Inlet: 1999) and submarine (e.g., Kitimat Arm: 1975) failures, larger continental slope failures off the west coast of Vancouver Island, and those associated with the deformation front of the Cascadia subduction zone, with the potential to increase the size of megathrust-generated tsunamis, are the most pertinent potential sources of a landslide tsunami for the west coast of Canada (Clague et al. 2003;Goff et al. 2020;Leonard et al. 2014;Murty 1979). Potential far-field landslide tsunami sources include the Aleutian trench and volcanic flank collapses in Hawaii or other Pacific islands (e.g., Leonard et al. 2014).
The shores of B.C. have the highest tsunami hazard in Canada, with a significant population at risk, particularly along coastlines of the Strait of Georgia and Juan de Fuca Strait (which together with Puget Sound comprise the Salish Sea), as well as at the heads of inlets (e.g., Kitimat). This region also features a complex maritime environment as well as high onshore terrain, local steep slopes, and tidal variation, all of which contribute to a high risk of damage from subaerial and submarine landslides and related tsunamis. The damage risk increases if paired with high tide and/or storm surges.
The hazard of landslide tsunamis on the west coast of Canada remains poorly understood; the historical and geological record of slope failures is incomplete, as much of the seafloor remains to be mapped or studied in detail (Lintern et al. 2020). Only a few landslide tsunami modeling studies have been completed for the region, including the 1975 submarine failure in Kitimat Arm (Skvortsov and Bornhold 2007), the sixteenth-century subaerial failure in Knight Inlet , and potential submarine failures of the Fraser delta front (Dunbar and Harper 1993;Rabinovich et al. 2003), and offshore Texada Island (Rabinovich et al. 2003).
This study focusses on the potential for a large tsunamigenic subaerial failure off Orcas Island in the Strait of Georgia (Fig. 1). This work forms part of a larger project to assess coastal flooding hazards in Boundary Bay, near the US/Canada border, from a variety of sources. Although no tsunamigenic landslides have been documented in the short written historical record (∼< 200 yr), recent landslide activity is indicated by a rockfall debris apron as well as two distinct landslide deposits on the seafloor adjacent to northeast (NE) Orcas Island ( Fig. 1; Greene et al. 2018).
Both slides initiated subaerially on the steep NE side of the island, and one of the slides has a 450-m runout (Greene et al. 2018). A likely trigger for future slides here is strong ground shaking during large earthquakes on the nearby Skipjack Island fault zone, which has been active in the Holocene and may host earthquakes up to magnitude 7 or higher (Caston 2021;Greene et al. 2018). A sudden collapse along the steep NE coast of Orcas Island is expected to cause tsunami waves with the potential to strike nearby coastal areas of the Gulf and San Juan islands and the adjacent mainland-modeling is needed to assess this potential.
We carry out numerical simulations of a potential subaerial landslide on NE Orcas Island and the resultant tsunami waves to simulate their potential impacts on nearby coastlines in the southern Strait of Georgia near the US/Canada border (Fig. 1). Using ArcGIS 10.7 software, we first compute the physical properties of the subaerial landslide (SAL), such as geometry (length, width), center of mass, and volume (Sect. 2.1). The kinematic characteristics of the landslide, such as initial acceleration and final velocity, are calculated using the law of motion in Sect. 2.2. Simulation of the subaerial landslide-generated tsunami (SALGT) is then performed using NHWAVE software (Ma et al. 2012) (Sects. 2.3 and 2.4), and the results of the numerical simulations are interpolated to the FUNWAVE-TVD (Total Variation Diminishing) variant of the fully nonlinear Boussinesq wave model (FTVD) (Kirby et al. 2013;Shi et al. 2012). This model simulates tsunami wave propagation toward the shoreline (Sects. 2.3 and 2.4), enabling us to determine inundation, runup heights, and current velocities that are described in Sect. 3. Fig. 1 a Map of the study area, highlighting the locations of the landslide source (Orcas Island) and Boundary Bay as the focus for assessing tsunami impacts; b evidence of subaerial landslide failures near Orcas Island (modified after Greene 2019); and c high-resolution multibeam echosounder (MBES) bathymetric and land LiDAR images indicating the area of the Skipjack Island fault zone and evidence of past landslides along northeastern Orcas Island (modified after Greene et al. 2018) 1 3 2 Subaerial landslide-generated tsunami (SALGT) modeling

Location and volume of the potential landslide
The boundaries, form, and volume of a landslide must be approximated for accurate numerical modeling. Subaerial slides have a greater tsunamigenic potential than submarine slides of equal volume and shape (Fornaciai et al. 2019). The most significant factor affecting the evolution of a slide-generated tsunami and associated wave heights is the slide volume (Chen et al. 2014;Løvholt et al. 2020).
In our modeling, we consider the largest likely tsunamigenic subaerial landslide on NE Orcas Island. Given the occurrence of two large landslide deposits on the adjacent seafloor and a debris apron that extends between the deposits, we postulate that the entire steep slope of NE Orcas Island between the two large deposits is susceptible to a future failure. This failure is considered to be a worst-case scenario. To pinpoint the position of the probable landslide, geological mapping data were merged with digital elevation models (DEMs) obtained from multibeam bathymetric and topographic survey data. The potential failure area is shown in Fig. 2; its cross-slope width and downslope length are estimated to be ∼ 5840 m and 1080 m, respectively.
The volume of the potential landslide was calculated by approximating a post-failure surface that is based on an interpolation between the existing large landslide deposits and analyzing the volume difference between the pre-and post-failure DEMs. The volume was computed using the ESRI ArcGIS 10.7 GIS software's "Spatial Analyst Cut Fill" technique to calculate projected net-loss volumes. This method generates a map based on two input raster surfaces of the potential slide area, one before and one after the slide occurrence; the seafloor map is characterized by the difference between the pre-slide bathymetry and a virtual plane (the approximated post-failure surface) (e.g., Verbovšek et al. 2017). The difference between the pre-and post-failure surfaces on the seafloor is then subtracted from a digitized polygon representing the onshore potential failure area in order to produce the post-failure subaerial surface. We used a global polynomial trend interpolation set with a polynomial order of three (P3) to create the virtual plane surfaces. The pre-and post-failure surfaces were then run through the "Cut Fill" processing tool to calculate the volume of the potential slide (Verbovšek et al. 2017). The tool compares the heights of raster values in each cell. The volume difference is given as positive if the material was removed (blue in Fig. 2c) or negative if the material was added (filled) (red in Fig. 2c). This estimated volume approach does not take into account the other more complex strategy of anticipating volumes for un-failed slides, e.g., using data from core profiling or seismic surveys, as no such data are available. The volume of the potential worst-case landslide is estimated at 0.17 km 3 (Table 1).

Kinematic and geometric parameters of the subaerial landslide
As previously noted, landslide volume is the most important factor determining tsunamigenesis, but other variables such as landslide frontal area and landslide dynamics including velocity, Froude number, and acceleration also play a significant role (Løvholt et al. 2020).
Sliding and slumping are the two most common mechanisms for subaerial and submarine mass failures (MF), as rigid or deforming translations down a flat slope . The MF center of mass motion is determined by geometric, hydrodynamic, and geotechnical characteristics, depending on whether the MF is rigid or deforming .
Using the estimated geometry of the Orcas Island SAL, we numerically modeled the tsunami generation and propagation. No mass (material) was extracted from the bathymetric data before the simulations. We assume a worst-case failure scenario consisting of a slump, a rigid and impulsive subaerial landslide. Rigid landslides are more tsunamigenic and create more inundation than deformable landslides (Schambach et al. 2019). In addition, a rigid landslide can be characterized in numerical simulations as a mound with a "Quasi-Gaussian" cross-sectional shape of elevation ( , , t) , an elliptical footprint of downslope length "b" and cross-slope width "w", and a maximum thickness "T" (Grilli et al. 2015).

3
The geometric parameters and beginning position determined for the modeled landslide are given in Table 1 and Fig. 3. We can only make estimates about the probable landslide thickness due to a lack of geophysical and morphological knowledge about the landslide site. The maximum thickness, T = 75.7 m, is determined using the volume produced from ArcGIS 10.7 software, the elliptical footprint of the landslide (Table 1), and the V s = 0.351 bwT equation (Grilli et al. 2015). These dimensions provide a T/b value of 0.07, which lies at approximately the lower end of the T/b validity range of [0.1, 0.15] determined by Watts et al. (2005) based on 12 numerical simulations of underwater slumps.
The coordinate transformation to the local MF slope-parallel coordinate system ( , ) is specified based on the MF starting center of mass position (Table 1), in global axes (x, y) and azimuthal direction of motion, (as defined in Schambach et al. 2020). Grilli and Watts (2005) addressed 3D problems in the vertical plane for slices of width "w" since it is assumed that mass failures are homogeneous in width. Further assumptions include the absence of isolated high Reynolds number flows near moving MFs, which allows the fluid to be assumed to be inviscid and irrotational .
The kinematic properties of the rigid subaerial slump are calculated using the analytical rules given by Watts (1999, 2005) and Watts et al. (2005). The kinematics for the rigid body motion are identical to its center of mass motion, S c (t) (provided by Eq. (1) below), in azimuthal direction = 32 • on the surface of uniform slope s = 34 • . (Grilli et al. 2015).
Based on prior research (Locat et al. 2009), it can be shown that the maximum displacement (runout), S f , is rather small in its leading direction of movement down the slope ( S f < b ) throughout its tsunamigenic movement, spanning an indeterminate duration of activity, t f (Grilli et al. 2015).
U max = S 0 ∕t 0 is the maximum slump velocity (Schambach et al. 2020). Similarly, the slump acceleration A(t) is calculated using the second derivative of displacement based on Eq. (1), where A 0 is the starting acceleration (Grilli et al. 2015). U max , A 0 , and other kinematic parameters estimated for the Orcas Island SAL are given in Table 2.
Although landslide velocity is important for tsunamis generated by impulsive slides, such as fast-moving SALs, it is well recognized that landslide acceleration dictates the initial tsunami elevation generated by translational landslides (Løvholt et al. 2015). In other words, the starting acceleration, which impacts initial elevation in the event of a lengthy run-out, is the most important kinematic parameter for a tsunamigenic landslide (Løvholt et al. 2015). SALs, unlike totally submerged landslides, frequently entail significant Froude numbers and nonlinearity (Løvholt et al. 2020). The landslide Froude number also serves as the most important kinematic property for determining the height of the resultant tsunami for SALs and slumps where the runout distance is sufficiently short in relation to the slide length (Løvholt et al. 2015). As a result, the landslide moves impulsively, producing shorter and higher waves that are subject to increased frequency dispersion (Løvholt et al. 2015).
Due to the low velocity and modest amplitude of motion, little hydrodynamic drag is anticipated for rigid slumps (Grilli et al. 2015). Thus, inertia comprises both the MF mass M s = s V s , with s representing the slide bulk density, and the specific density = s ∕ w , with w denoting the water density, as well as an added mass ΔM s = C M w V s , defined by way of an added mass coefficient C M . For rigid slumps failing in a circular motion of high radius R and modest angular displacement ΔΦ , we can also assume constant basal friction. By formulating a dynamic balance between weight, buoyancy, inertia, and basal friction forces (Grilli et al. 2015) developed an analytical rule of motion for the center of mass: where "g" represents earth's gravitational acceleration. The radius of the slump motion is calculated using the semi-empirical formula provided by Watts et al. (2005) as a function of slump downslope length and maximum thickness.
We may deduce from Eq. (1) to (2) that the slump law of motion is a function of runout distance S f , and the slump time of motion t f is a function of = s ∕ w , the specific density (Grilli et al. 2015).
In our study, we assumed a slide material density of s = 1790 kg∕m 3 , an average sea water density of w = 1025 kg∕m 3 (so, = 1.75 ), a hydrodynamic added mass coefficient of C M = 1 , and an angular displacement of ΔΦ = 0.38 radians (or 21 • ) based on Watts et al. (2005).
Using Eq.
(2), a radius of nearly circular motion R, characteristic time t 0 , and characteristic distance of motion S 0 , are estimated. Then, the maximum displacement ( S f ) and runout time ( t f ) can be derived from Eq. (1), giving S f = 732 m and t f = 84 s in azimuthal direction = 32 • (Table 2).
Parameters computed in Tables 1 and 2 are used in the slump law of motion in addition to slump geometry to assign the bottom boundary condition in NHWAVE, described in Sect. 2.4 (Grilli et al. 2015;Schambach et al. 2019Schambach et al. , 2020.

Computational grids, and bathymetry/topography data
Based on guidelines from the US National Tsunami Hazard Mitigation Program (NTHMP) (Dunbar and Weaver 2015), the resolution in coastal regions for evaluating coastal runup should be 10 m or better. Use of a coarser resolution grid along the coast cannot properly predict tsunami coastal inundation and runup; e.g., coarse resolution may create obstructions that prevent waves from entering a port and limit the maximum wave height and inundation (Nemati et al. 2019).
In this study, we integrate and interpolate the best available bathymetry and topography data into five horizontal Cartesian grids of different resolutions ( Fig. 4; Table 3). We used 10 m topo-bathymetry data from Natural Resources Canada (2021) to create a seamless topographic/bathymetric DEM for the grid G0 region with 160-m resolution. For the other grids, G1, G2a, G2b, and G3N, we utilized data from a cross-border 3-m resolution DEM compiled by various government agencies (Amante et al. 2020;Carignan et al. 2019). The resulting bathymetry and grid areas are shown in Fig. 4. G0, G1, G2a, and G2b nested grids, with resolutions of 160 m, 40 m, 10 m, and 10 m, respectively, are utilized in the FTVD simulations ( Fig. 4; Table 3). The grid G0 covers the largest area, including northern Washington State, southern Vancouver Island, eastern Juan de Fuca Strait, southern Strait of Georgia, and the adjacent B.C. mainland; G1 encompasses the southernmost portion of the Strait of Georgia and adjacent parts of the mainland and Vancouver Island. The two 10-m resolution grids, G2a and G2b, are nested in order to more rigorously simulate tsunami inundation and runup in two areas of particular interest: (1) the coastal parts of Orcas Island, and (2) the Boundary Bay region. G3N, also at 10-m resolution, covers the eastern sector of Orcas Island and is used to model landslide generation in the NHWAVE model, with five evenly spaced boundary fitted layers in the vertical (from the sea surface to the seabed).
The effect of ocean tides on inundation and runup is another important component in tsunami modeling; a higher tide at the time of the landslide leads to greater inundation for the same tsunami source. We used a tide level of Higher High Water Mean Tide (HHWMT), the average of all Higher High Water projections at a site over the previous 19 years; this is equivalent to Mean Higher High Water (MHHW) in the USA, typically used in tsunami modeling (e.g., Suleimani et al. 2005). The reference level in the topobathymetry dataset (grid G0) is chart datum (CD), roughly comparable to the lowest astronomical tide level (Vancouver Fraser Port Authority 2020); thus, the dataset must be corrected to HHWMT.
The reference level of topo-bathymetry data for all the other grids (G1, G2a, G2b, and G3N) is based on the Canadian Geodetic Vertical Datum 2013 (CGVD2013), which is deemed comparable to Mean Water Level (MWL), but with an offset of − 0.2 m. The next step in the modeling was to add the Higher High Water Mean Tide (HHWMT) level to CD or CGVD2013 levels in order to incorporate the static tide influence into the topo-bathymetry data. In this case, we used the HHWMT from the Canadian tide and current tables (Canadian Hydrographic Service 2022), which is based on CD, for the Fulford Harbour and Point Atkinson stations, which are the closest stations to the study area (see Fig. 4 for locations). Fulford Harbour, on the southeast coast of Salt Spring Island, B.C., has an HHWMT value of 3.3 m, whereas Point Atkinson, in West Vancouver, has an HHWMT value of 4.5 m. Grid G0 bathymetry is updated with the average HHWMT of these stations (3.9 m). The other grids are based on MWL, so must be corrected for the difference between MWL and HHWMT. This difference, the tidal amplitude, is 1 m at Fulford Harbour and 1.4 m at Point Atkinson; we use the average of 1.2 m to correct the reference level of the topo-bathymetry data. This result is also consistent with prior research conducted by Fine and Thomson (2020).

Numerical models and tsunami modeling methodology
Tsunami waves generated by landslides are more challenging to model than those generated by earthquake displacements (Løvholt et al. 2015). Landslide collapse at the water body boundaries causes impulsive waves that spread outward and migrate along shorelines (Romano et al. 2019). The relative wave amplitude of a SALGT depends on both the Froude number (F) and the slide thickness, as confirmed in numerous studies (e.g., Fritz et al. 2004;Heller and Hanger 2010;Fornaciai et al. 2019;Sarlin et al. 2021).
The tsunamigenic potential of a slide is determined by both the Froude number (F) and the slide thickness (Fornaciai et al. 2019). In order to determine how effectively landslides generate waves, we must calculate the Froude number, which is the ratio of landslide speed to wave speed (linear shallow water wave celerity = √ gh where h is the depth of the water). Wave generation is most efficient if the landslide and tsunami waves travel at the same velocity, whereby F = 1 . As the Froude number decreases below unity, wave generation becomes less efficient (Fornaciai et al. 2019). In the case of the Orcas Island landslide, the maximum velocity of the landslide as it moves downward to the sea is around 13.64 m s −1 , and the water depth adjacent to the failure area is 75 m. Consequently, the Froude number is about 0.5, indicating moderately effective wave generation due to the landslide.
To model the Orcas Island landslide and subsequent tsunami, we use a coupled modeling approach, similar to that used in recent studies on landslide processes (e.g., Grilli et al. 2013Grilli et al. , 2015Schambach et al. 2019Schambach et al. , 2020Schambach et al. , 2021. Based on the dimensions determined in Sects. 2.1 and 2.2, a simplified landslide form and plane geometry were computed on grid G3N (Fig. 3).
The landslide failure movement and resulting tsunami production are modeled using the three-dimensional (3D) non-hydrostatic numerical wave model NHWAVE ). In the horizontal direction, the model employs a Cartesian system (Δx, Δy) , and in the vertical direction it uses a -layer boundary fitted grid (Schambach et al. 2020). With fully nonlinear free-surface boundary conditions, NHWAVE solves the inviscid Euler equations (viscous and turbulent effects can be incorporated if needed) . In most cases, because tsunami waves are approximately 2D, only a few -layers are needed (Schambach et al. 2020). Adding more -layers better quantifies wave dispersion, especially where the depth to wavelength ratio is relatively high. For coastal landslide tsunami simulations, five -layers appear to suffice (Grilli et al. 2015;Schambach et al. 2019Schambach et al. , 2021. The use of three -layers in NHWAVE results in the same dispersive features as a 2D Boussinesq model (Schambach et al. 2020).
In the near field, five boundary-fitted vertical layers in the G3N grid are taken to estimate the lateral collapse landslide and the resulting tsunami production in the NHWAVE setup. Using NHWAVE to simulate wave generation helps us to improve the assessment of dispersive effects and horizontal velocities. In landslide tsunami generation and propagation models, the availability of a frequency dispersion option has a substantial impact. Dispersion governs phase speed, wave-wave interactions during propagation, and tsunami impact (Schambach et al. 2021). Undular bores (also known as dispersive shock waves) are thought to be caused by dispersion across the crest and trough of prolonged shoaling tsunami waves near coasts, boosting coastal impact (Schambach et al. 2021;Madsen et al. 2008).
The modeled landslide occurs on the steep NE coastline of Orcas Island. When tsunami waves occur in shallow water, the relationship between the waves and the sloping seafloor is a crucial element to consider (Romano et al. 2019). The waves can be refracted or trapped (Romano et al. 2019). Therefore, the simulation begins on grid G3N (10-m resolution) with the most recent version of NHWAVE (v.3.0), which can evaluate vertical acceleration impacts in the slide layer (Zhang et al. 2021a, b); this NHWAVE advancement is especially important on steep slopes Schambach et al. 2021).
Once the tsunami waves are generated, the waves are input to nested grids of increasing resolution toward the shore using the 2D fully nonlinear and dispersive Boussinesq wave model FTVD . The model simulates propagation from the near field to the surrounding coastlines in the far field.
Both directly forced waves and freely propagating waves occur during the landslide. The waves radiating away from the slide will encounter modified depths associated with the slide volume. It is these waves that need to take into account the changes in seafloor depth due to the slide. Unlike FTVD, NHWAVE can simulate moving bottom-generated waves in three dimensions at each time step, allowing velocities to be captured with more precision, as they change with depth during wave formation. As a result, landslide tsunami generation may be precisely simulated. On the other hand, coastal wave changeovers, such as wave breaking and dissipation, as well as incorporation of a shifting coastline, can be evaluated more precisely in FTVD. Furthermore, the use of a 2D grid in FTVD means that the computational time is a factor of three less than the required time for an identical horizontal resolution for a NHWAVE simulation. Owing to its spherical implementation, FTVD may be used for far-field tsunami propagation, whereas NHWAVE can only examine a Cartesian horizontal grid, limiting its applicability to near-field tsunami generation. These are the justifications for using linked modeling approaches in our research (more detail in Grilli et al. 2015). NHWAVE and FTVD are both open source (available on Github) and parallelized with Message Passing Interface (MPI) to operate effectively on big computer clusters (Schambach et al. 2020). According to NTHMP standards, the model has been fully verified based on laboratory and field benchmarks for both coseismic and rigid landslide tsunami generation and propagation Tehranirad et al. 2012).
A similar coupled modeling strategy was previously used to simulate the coastal effect of transoceanic tsunamis along the US east coast, such as those resulting from a potential collapse of the Cumbre Vieja volcano in the Canary Islands (Abadie et al. 2012), for which the SAL tsunami source was computed using the multi-fluid 3D Navier-Stokes solver THETIS (Abadie et al. 2010). Grilli et al. (2013) used an equivalent NHWAVE/FTVD coupling approach to model the Tohoku 2011 tsunami, where the seismic source was located as a time-and space-varying bottom boundary condition.
The Orcas Island landslide-generated tsunami was modeled using the initial landslide geometry and net-loss volume, kinematic parameters, and 10-m resolution bathymetry data. The modeled landslide stopped moving after runout time t f = 85 s, and modeling of the tsunami waves was continued in NHWAVE for an additional 15 s. There is no need to continue modeling in NHWAVE for a longer time, since landslides are not tsunamigenic after they have stopped moving. An animation of the NHWAVE model simulation (Vide-oNo01.mp4) is provided in Online Resource 2.
After 100 s, tsunami simulation is continued using FTVD, with initial conditions based on the NHWAVE simulation results. Surface elevation, horizontal velocity, and afterfailure bathymetry are interpolated in grid G0, G1, and G2a. Tsunami propagation is then simulated from the source location to surrounding coastlines in a succession of layered and nested grids, with increased resolution in the areas of particular interest, as outlined in Sect. 2.3.

Absorbing (sponge) layers define open boundary conditions in both models.
To prevent reflection at open boundaries, sponge layers are allocated, being 5 km wide at the west and 1 km at the south and north boundaries of grid G0. These widths appear adequate, as no reflection occurs at the boundaries in our results.
FTVD can compute inundation and runup with high precision at coastal locations by taking into account the impacts of bottom friction, dissipation, wave breaking, and the moving shoreline algorithm in the model . Where water depths are shallow, such as between Orcas Island and Boundary Bay, landslide-generated tsunami waves may not be widely dispersed in the beginning. The dispersive impact is also reduced if the landslide footprint is large in relation to the depth Schambach et al. 2019Schambach et al. , 2021. The footprint of the landslide in grids G2a and G2b is substantially larger (5.8 km width and 1.08 km length) than the maximum depth, which is roughly 180 m inside these grids. As a result, minimal wave dispersion occurs in those grid simulations.
In deeper water, as waves migrate away from shallow water locations, dispersion becomes more important. Significant dispersion occurs for the deeper-water grids G0 and G1.
We use a number of grid boundary stations (details not provided here) and tide gauge stations (Table 4; Figs. 5,6) to ensure data coupling between the grids. The boundary conditions for finer grids are based on the grid data at a time zero, i.e., the time series of boundary station elevation and horizontal velocity computed in grid G0 are given as the boundary condition for In the NHWAVE and FTVD model simulations, bottom friction associated with coarse material is assumed ( c d = 0.0025).

Table 4
Locations and modeled tsunami wave characteristics at numbered stations for which surface elevation time series are shown in Fig. 6 Station no.  1 3

Tsunami modeling results
The modeled tsunami waves generated by the potential Orcas Island landslide are shown and described in this section. Water surface elevations are shown as time series for the 12 selected stations in Fig. 6 (with plots for a total of 34 stations provided in Online Resource 1), and on maps in Figs. 7 and 8. Maps of current velocities are shown in Figs. 9 and 10 and the maximum tsunami runup within the high-resolution grid areas in Fig. 11. Table 4 contains characteristics of the modeled tsunami waves at the 12 selected stations. The 3D movement of the Orcas Island SAL is simulated in the NHWAVE model; the slide material flows downslope to the ocean, where the water depth is 75 m. During the northeastward propagation of the tsunami waves, the water depth remains uniform over a 1.7 km distance (Fig. 5). The depth then shallows to 40 m for about 1.4 km near Clark Island, before returning to a consistent depth of 100 m for another 2.16 km (to Lummi Island).
Surface elevation estimates from NHWAVE and FTVD up to t = 100 s show good agreement, especially for the crest-to-trough heights of the leading waves. This is an indication that the two models are well coupled. In grid G2a, immediately adjacent to the simulated subaerial landslide site, a huge initial tsunami wave is generated, with a wave crest amplitude of ∼ 20 m (Figs. 6-station 07, 7-b). The initial crest is followed by a ∼ 12-m trough. Tsunami wave reflection and shoaling are  Maximum modeled runup (color scales in meters) in map view for grids G2a (a) and G2b (b), and plotted against longitude for grids G2a (c) and G2b (d) primarily regulated by nearshore bathymetry, and the details of the tsunami source fade as the waves move from the source region . Hence, in the near-source region, the tsunami waves have very small wavelengths, with a period of approximately 400 s, and dissipate to near-zero amplitudes after 40 min due to a near-uniform depth, relatively high bottom friction, and marked dissipation along the shoreline (see Fig. 6-Stations 07, 08). As was previously indicated, a short runout distance with respect to the slide length may be the cause of the shorter wave duration, higher waves, and fast wave dissipation. Although the majority of wave energy propagates away from Orcas Island with the leading wave, some tsunami waves are re-directed by bathymetry and surrounding islands to eventually run up on shoreline areas around Orcas Island, away from the source location.
In the Eastsound and Westsound regions (see Fig. 5 for locations), the leading wave is 0.4 m and subsequent wave surface elevations increase to 1.24 m (detailed time series in the Online Resource 1, station 23); this occurs due to shoaling and wave buildup inside these semi-enclosed shallow areas that can occur due to the potential resonance of trapped waves (Nemati et al. 2019). Nearshore amplification effects are common and are often connected to bathymetric control (producing wave convergence/divergence) and coastline form (i.e., causing potential resonance) (Nemati et al. 2019).
In Grid G2b, the landslide-generated waves arrive at the southern section of the Fraser River Delta, south of Point Roberts, 13 min (780 s) after the landslide collapse, with roughly 0.3 m wave amplitude (Fig. 7). The greatest wave amplitude near the shoreline of Boundary Bay (Fig. 6-stations 09, 10) is around 1.2 m, with a 15-min arrival time. The highest wave elevation between White Rock and Semiahmoo Bay is around 1.04 m, with arrival times of about 25 min (more detail in Online Resource 1). The tsunami reaches Birch Bay (Figs. 5, 6-station 11, 7) within 20 min; the first two waves have amplitudes of 1-1.1 m. The region near Lummi Bay (e.g., Figs. 5, 6-station 12, 7) has the largest maximum amplitudes in grid G2b, of about 1.2 m, with a 12-min arrival time. In semi-enclosed locations, such as harbors or shorelines where the geometry is complicated, the reflected waves can undergo many oscillations. For example, at stations 09 and 10, the leading waves (max 0.45 m) are smaller than the crests of later waves (max 1.14 m); at these locations, the maximum wave height occurs 27 min and 10.3 min, respectively, after the leading wave arrival (Fig. 6). The explanation for this might be the location of the Serpentine River on the NE side of Grid G2b, where waves propagate to a shallower location and are not slowed by dry regions near Boundary Bay. Also, tsunami travel time to shore is governed by wave celerity, which is determined by water depth, and landslide-generated short waves are controlled by wave frequency (e.g., Grilli et al. 2019). The large leading waves at nearshore stations in grid G2b are accompanied by oscillations with a short average period of 400 s. This suggests that dispersion has less influence in shallow water zones. The dissipation begins after 100 min at all nearshore stations, and the wave heights gradually decrease.
In grid G1 (Fig. 5), the first wave arrives near Victoria, after 20 min, with the largest amplitude of 0.3 m occurring 6 min later (detailed time series in Online Resource 1). South of Orcas Island, the maximum wave surface elevation is less than 0.5 m. At station 03, located NW of Orcas Island near Pender Island (Fig. 5), the leading tsunami wave arrives within 15 min with an amplitude of almost 0.5 m; wave amplitudes decrease over time (Fig. 6). Within grid G1, in areas removed from the source and from the direct path of wave propagation, maximum wave amplitudes do not exceed ∼ 0.5 m. In contrast, the leading wave at station 04 near Lummi Bay (Fig. 5) has an amplitude of ∼ 4.4 m, with a subsequent 1-m trough, followed by smaller waves (Fig. 6). This station is positioned in the direct path of the waves generated by the tsunami source, with no dry regions or islands to scatter or deflect the wave motion. The wave amplitude reduces as it approaches the shore due to increased bottom friction and energy loss via breaking.
The time series at stations 01 and 02, located in the northern (01) and southern (02) areas of grid G0, ∼ 105 km and ∼ 65 km from the source, provide examples of the minimal far-field impact of the tsunami (Figs. 5, 6). At both sites, the leading wave arrives after ∼ 40 min, and wave amplitudes reach 0.14 m.
The propagation of the tsunami waves is illustrated for the high-resolution grids as a time sequence, up to 100 min following the initial landslide, of surface elevations in Fig. 7 and current speeds in Fig. 9. Full animations are also provided in Online Resource 2. Over the first ∼ 100 s, large leading waves (20 m high) propagate in all directions, but primarily toward the NE. At time t = 80 s (near the t = 85 s landslide runout time) tsunami waves are travelling in several directions, resulting in a complex tsunami wave pattern. Interactions between waves, bathymetry, the shoreline, and neighboring islands cause multiple wave reflections. Figure 8 shows the maximum tsunami amplitudes from the Orcas Island simulation. The largest amplitudes generally occur close to the Orcas Island source, with amplitudes generally > 15 m. The maximum wave height is 19 m at 50 s after landslide initiation. In the direction of wave propagation, including the southern Boundary Bay region, particularly from Lummi Bay to Neptune Beach (Fig. 5), amplitudes exceed 4 m (max. 4.41 m at 9 min). The region is distinguished by its numerous small inlets and semi-enclosed embayments, where tsunami waves can be amplified and potentially cause more damage. A reduction in water depth that causes wave shoaling and/or resonance, as well as limiting area, can contribute to wave amplification.
The strongest currents, up to 10 m/s (e.g., 9 m/s at station 07), coincide with the largest wave amplitudes in the near-source region . Strong modeled currents occur in coastal areas near the source region, as well as further to the NE within Boundary Bay, Semiahmoo Bay, and Birch Bay. Simulated currents are less than 4 m/s elsewhere in the region and are relatively weak in more remote areas of the model domain. Specifically, tsunami-induced currents throughout Boundary Bay have typical speeds of less than 1-1.5 m/s and do not exceed 3 m/s. Figure 10 depicts the maximum modeled current velocities. Fast current velocities (e.g., higher than 1.5 m/s) are forecast over the whole east, north, and south coast of Orcas Island in grid G2a (Fig. 10c). However, the results in the south part (Shaw Island) may be artefacts because they are at the grid's edge. The highest current velocities from Lummi Bay to Boundary Bay can be seen in Fig. 10d. Due to shallower depths and the narrow constrictions between the many islands in the direction of wave propagation, particularly fast, hazardous currents are predicted near Point Doughty, Buckhorn, and Clark Island in grid G2a, and Neptune Beach, Birch Bay, Semiahmoo Bay, White Rock, Boundary Bay, Point Roberts, Tsawwassen, and Richmond in grid G2b.
Maximum tsunami runup is shown in Fig. 11 for the two high-resolution areas of interest. Starting from the north and moving clockwise (toward Buckhorn and Doe Bay), the runup is initially high in the near-source region, ∼ 30 m, then lowers to less than 5 m before dramatically increasing to a local maximum of 30 m on Shaw Island, surprising given its location far removed from the source area (Fig. 11a, b). This extreme runup might be caused by being located on the boundary of computing grid G2a, which could lead to some inaccuracy and artificial amplified wave height, causing a spurious runup prediction in this area. Due to refraction surrounding the landslide, runups on Point Doughty, which is located approximately perpendicular to the primary collapse direction, still reach 10-16 m. Smaller, but still significant, runup occurs throughout the region surrounding Orcas Island. Further afield to the NE (Fig. 11c, d), runup reaches 7.5 m at Neptune Beach near Lummi Bay. Both initial and reflected waves cause significant runup (> 1.5 m) along much of the shoreline between Point Roberts and Lummi Bay.

Discussion
In this study, we simulated tsunami waves triggered by a potential subaerial rigid landslide on NE Orcas Island, using NHWAVE ) and FTVD (Kirby et al. 2013;Shi et al. 2012) software.
Landslide extent and volume were determined based on a worst-case scenario: largescale failure of NE Orcas Island for the full length between previous landslide deposits. Small differences in volume or rheology are expected to have little effect on a tsunami coastal impact projection . The sensitivity of the numerical results to model formulation, seabed topography, bottom friction, and other physical factors, as well as to tsunami wave resonances and attenuation time scales, requires more investigation.
Grids G0 and G1, which include areas of relatively deep water, were used to run simulations with dispersive effects. Initial testing showed that it was not necessary to include dispersion for simulations in the shallower regions of the G2a and G2b grids. All simulations in FTVD were performed out to t = 7200 s (120 min), which was long enough to capture the maximum runup throughout the Orcas Island and Boundary Bay regions.
Tsunamis generated by submarine landslides begin with a trough (negative wave), whereas subaerial landslide-generated tsunamis (SALGTs) commence with a crest (positive wave) followed, as in Fig. 6, by a trough of similar or greater magnitude (Fornaciai et al. 2019). According to Løvholt et al. (2015), the leading wave tsunami amplitudes caused by SALs decay more slowly than those caused by totally submerged landslides. The far-field propagation of tsunamis caused by rockslides in fjords has comparatively unique properties compared to offshore occurrences, due to the multiple reflections that make the propagation less efficient (Løvholt et al. 2015). Except for the near-source region offshore Orcas Island (stations 07 to 08), these characteristics of SALG waves are plainly visible at stations within our study area (Fig. 6). Dispersion and geometric spreading may be the two primary explanations for the relatively rapid decrease in surface elevation near the tsunami source for stations 07 and 08 during wave propagation. As the waves travel to the northeast from the source, these variables, together with the transfer of tsunami energy into channels, result in rapid attenuation of the waves. These factors can also have an influence on runup.
The effect of waves on the shore is not only directly connected to the source from which they originate, but also highly influenced by topography and bathymetry. The distance from the tsunami source does not have a linear connection with wave height or runup. Bathymetry is critical since it may lead to a significant increase in wave height in a given region. Topography, which is equally crucial, defines the inundation zone.
The highest runup is captured by the finer grids G2a,b rather than the coarser grids G0 and G1 (more detail in Online Resource 3), justifying the use of nested grids in FTVD. As shown in Fig. 11, the expected tsunami impact is considerable along the east coast of Orcas Island, near the source, and for much of the mainland coast between White Rock, B.C. and Lummi Bay, Washington, an area that is positioned directly in the tsunami propagation path, without intervening islands to cause damping. Elsewhere, wave heights (and runup) are damped by islands or diminished with distance due to dispersion and geometric spreading.
We can compare our numerical modeling results with the maximum near-field wave amplitude predicted by experiment-based equations. Mulligan and Take (2017) first examined the condition of the bulk fluid upon landslide impact using a one-dimensional technique. The hydrostatic pressure gradient formed by a near-instantaneous vertical waterlevel displacement over a short timescale Δt e = t(g∕h) and length scale L e expresses the momentum imparted to the water in this situation. The initial pressure gradient is calculated from the difference in water surface elevation (x, t) between two locations in the landslide impact zone and the zone unaffected by the slide, which has a maximum positive wave amplitude above the still water depth h. The analysis yields where a q denotes the maximum near-field wave amplitude derived from a quadratic equation, h is the still water depth, s is the sediment density, T is landslide thickness, is slope angle, u max is maximum velocity parallel to the slope, the fluid density, and g is gravitational acceleration (Mulligan et al. 2018). In the case of the Orcas Island landslide, the near-field water depth is 75 m, and the other parameter values are as follows: s = 1750 kg∕m 3 , u max = 13.64 m∕s, = 34 • , L e = 732 m, Δt e = 85 (calculated based on period of the first crest in station 07 as 10 s) which gives a maximum near-field amplitude of about 15.8 m. This is slightly smaller than the maximum amplitude of ∼ 19 m in our model results (e.g., Fig. 6-station 07).

Conclusion
In this study, we used the dispersive numerical tsunami models NHWAVE  and FTVD (Kirby et al. 2013;Shi et al. 2012) to simulate tsunami waves generated by a large potential subaerial landslide on the NE coast of Orcas Island.
Modeling of the subaerial landslide (SAL) necessitated the use of high (3 m)-resolution bathymetry data (cf., Amante et al. 2020;Carignan et al. 2019). The simulated SAL travelled around 732 m before coming to rest in water depths of around 75 m. Based on an assumed friction drag coefficient of 0.0025 between the slide and the underlying seafloor, the modeled slide moved downslope with a peak velocity of 13.64 m/s before coming to rest after a duration of about 85 s. According to our numerical simulations, the SAL would generate tsunami waves with peak amplitudes of 15-20 m, current speeds of up to 9 m/s, and wave periods of roughly 6.6-8 min.
The modeled tsunami waves experience many reflections and a significant degree of dispersion due to their relatively short wavelengths and the complicated coastline and bathymetry in the Strait of Georgia. The most powerful waves and currents will occur at the shorelines adjacent to the failure on Orcas Island and opposite to it, near Lummi Bay. Wave amplitudes then rapidly attenuate with distance from the source due to the dispersion effect, geometric spreading, and transfer of tsunami energy in semi-enclosed areas. The wave amplitude near the tsunami source is about 17-19 m and the current speed and (3) a q = √ h 2 + 2 s Tu max cos L e gΔt e − h the maximum runup are roughly 9 m/s and 30 m, respectively. The simulated waves reach Lummi Bay (station 04, 15.5 km from the source region) in 8.5 min, with a propagation speed of approximately 31 m/s, maximum amplitude of ∼ 4 m, runup of ∼ 7.5 m, and current speeds of 2.17 m/s. Modeled tsunami waves reach the entrance of Boundary Bay in less than 12 min and arrive at Semiahmoo Bay in less than 15 min, with wave amplitudes of up to 1 m, accompanying currents of up to 1.5 m/s, and maximum runup of 4.5 m. Studies of past tsunami events (e.g., Melgar 2021; Whitmore et al. 2008;Williamson et al. 2020) demonstrate that coastal damage is possible even for tsunami amplitudes as low as 0.5 m, with more severe damage and inundation occurring for amplitudes or runup exceeding ∼ 1.5 m. Runup maps (Fig. 11) reveal that areas closest to the source area, including Point Doughty, Buckhorn, and Clark Island, would be the most impacted (runup > 12 m), while other locations on and near Orcas Island, including Eastsound, Westsound, and Waldron Island, would be less impacted with runup less than 5 m. The modeling shows that Clark Island, which hosts a marine state park and campground, could be completely overwhelmed by waves within ∼ 2.5 min of the landslide. High runups (up to ∼ 7.5 m) and strong currents (up to 4 m/s) would affect Semiahmoo Bay, Birch Bay, and Neptune Beach, which are located at the southern end of the Boundary Bay region. In contrast, Point Roberts, Boundary Bay, and White Rock would experience lesser wave effects with maximum runup of up to 3.5 m.