The Lituya Bay landslide-generated mega-tsunami – numerical simulation and sensitivity analysis

The 1958 Lituya Bay landslide-generated mega-tsunami is simulated using the Landslide-HySEA model, a recently developed finite-volume Savage–Hutter shallow water coupled numerical model. Two factors are crucial if the main objective of the numerical simulation is to reproduce the maximal run-up with an accurate simulation of the inundated area and a precise recreation of the known trimline of the 1958 mega-tsunami of Lituya Bay: first, the accurate reconstruction of the initial slide and then the choice of a suitable coupled landslide–fluid model able to reproduce how the energy released by the landslide is transmitted to the water and then propagated. Given the numerical model, the choice of parameters appears to be a point of major importance, which leads us to perform a sensitivity analysis. Based on public domain topo-bathymetric data, and on information extracted from the work of Miller (1960), an approximation of Gilbert Inlet topo-bathymetry was set up and used for the numerical simulation of the mega-event. Once optimal model parameters were set, comparisons with observational data were performed in order to validate the numerical results. In the present work, we demonstrate that a shallow water type of model is able to accurately reproduce such an extreme event as the Lituya Bay mega-tsunami. The resulting numerical simulation is one of the first successful attempts (if not the first) at numerically reproducing, in detail, the main features of this event in a realistic 3-D basin geometry, where no smoothing or other stabilizing factors in the bathymetric data are applied.

1 Introduction

Tsunamis are most often generated by bottom displacements due to earthquakes. However, landslides, either submarine or subaerial, can also trigger devastating tsunami waves. In addition, they are, on some occasions, extremely destructive as they form near the coast or in the same coastline if they are aerial. Sometimes landslides may generate so-called mega-tsunamis, which are characterized by localized extreme run-up heights (Lituya Bay, 1958, Miller, 1960; Fritz et al., 2009; 1934 Tafjord event, Norway, Jørstad, 1968; Harbitz et al., 1993, Taan Fjord, 17 October 2015, Bloom et al., 2016; Higman et al., 2018; or the very recent event of Sulawesi, 29 September 2018, among many others).

For seismic tsunami simulations, in general, the most critical phases are generation and arrival at a coast, including inundation. Propagation over deep basins can be modeled using the nonlinear shallow water (NLSW)  equations or more typically using a non-diffusive linear approximation. With landslide-generated tsunamis, however, matters become more complicated. The generation phase itself becomes critical and complex effects between the landslide and the water body must be taken into account. Most notably, the case of subaerial-landslide-generated tsunamis is where modeling and numerical implementation becomes most critical, owing to these events producing more complex flow configurations, larger vertical velocities and accelerations, cavitation phenomena, dissipation, dispersion, and complex coupled interaction between landslide and water flow. It is evident that shallow water models cannot take into account and reproduce all of these phenomena, particularly vertical velocities or cavitation. We do, however, demonstrate here that such models can, despite limitations, be useful for hazard assessment in which the main features (from a hazard assessment point of view) of these complex events, such as run-up and main leading wave, are reproduced. The overall aim is not to accurately reproduce the evolution of the displaced solid material or the dispersive nature of the trailing waves as the perturbation propagates, but, alternatively, to accurately reproduce the impact of tsunami waves to coasts in terms of run-up and flood area. Comparison of the numerical results with the observed trimline presented here is shown to support our statement that a fully coupled, vertically integrated shallow water Savage–Hutter model can, effectively and accurately, reproduce the run-up and coastal inundation resulting from aerial landslides generated in fjords and enclosed basins. This study further supports this assessment by comparing model results with observed data for a paradigmatic example of extreme run-up produced by an aerial landslide in an enclosed bay, a simulation that has not been successfully undertaken previously with more comprehensive numerical models.

Full 3-D numerical modeling of landslide-generated tsunamis (Horrillo et al., 2013) is uncommonly used for real-world scenarios due to the highly demanding computational resources required. Whereas common thought persists that NLSW equations are sufficient for simulation of ocean-wide tsunami propagation averaged models, the importance of frequency dispersion for modeling landslide-generated tsunamis lends preference to Boussinesq models. In no case, however, can any of these standard models describe the violent impact of subaerial landslides with flow separation and complex subsequent flow patterns and slide material evolution. We further claim that, in enclosed basins, the two main mechanisms that need to be well captured and accurately reproduced by a numerical model are first the transmission of energy from the slide material into the water body and, second, the coastal inundation by means of accurate wet–dry treatment. Confinement makes the role of dispersion minor relative to other effects.

In the case of tsunamigenic aerial landslides in fjords, bays, or any long and narrow water body, confinement and reflection (a process that also makes propagation and interaction more complex) are relatively more important considerations than dispersion, which becomes less important. This is particularly true for the leading wave (Løvholt et al., 2015) that, conversely, is mainly responsible for coastal impact. Lindstrøm et al. (2014) using a scaled laboratory setup showed that wave propagation along the fjord involved frequency dispersion but only to a moderate extent (Løvholt et al., 2015). It is in the far field that dispersive effects are proven to be important for a realistic description of tsunami impact (Løvholt et al., 2008; Montagna et al., 2011).

Despite this, and independent of the eventual confinement of the flow, many authors continue to claim the absolute need for dispersive, or even full Navier–Stokes, models (Abadie et al., 2009), when dealing with the simulation of landslide-generated tsunamis. Still other authors, not so strict in their assessment, claim that dispersive models represent a better alternative than NLSW models (Løvholt et al., 2015). The lack of dispersive model simulations in the literature of the Lituya Bay event plays against the argument of ruling out nondispersive coupled models, such as the one proposed in the present work. Very recently, in February 2017 during the “2017 US National Tsunami Hazard Mitigation Program (NTHMP) Tsunamigenic Landslide Model Benchmarking Workshop”, held on the Texas A&M Galveston, Texas, campus, participants agreed to recommend to the NTHMP Mapping and Modeling Subcommittee (MMS) the use of dispersive models for numerical simulation of landslide-generated tsunamis. For the case of enclosed basins, bays, and fjords, it was agreed that NLSW models remain a suitable tool.

Among all examples of subaerial-landslide-generated tsunamis, the Lituya Bay 1958 event occupies a paradigmatic place in the records, standing alone as the largest tsunami ever recorded and representing a scientific challenge of accurate numerical simulation. Based on generalized Froude similarity, Fritz et al. (2009) built a 2-D physical model of the Gilbert Inlet scaled at 1 : 675. A number of works have focused their efforts in trying to numerically reproduce Fritz et al. (2009) experiments (Mader and Gittings, 2002; Quecedo et al., 2004; Schwaiger and Higman, 2007; Weiss et al., 2009; Basu et al., 2010; Sánchez-Linares, 2011). Detailed numerical simulations of the real event in the whole of Lituya Bay with a precise reconstruction of the bottom bathymetry and surrounding topography are limited. As an example, Mader (1999) fails in reproducing the 524 m run-up and concludes that the amount of water displaced by a simple landslide at the head of the bay is insufficient to cause the observed tsunami wave. As far as we know, the present work represents the first successful attempt to realistically simulate and reproduce the 1958 Lituya Bay mega-tsunami in a realistic three-dimensional geometry with no smoothing in the geometry or initial conditions.

2 Background

At 06:16 UTC on 10 July 1958, a magnitude M_w 8.3 earthquake occurred along the Fairweather Fault (Alaska, USA). This quake triggered a landslide of approximately 30.6 Mm³ in Gilbert Inlet (Miller, 1960) that in turn produced the largest tsunami run-up ever recorded (Fritz et al., 2009). The epicenter of this quake was a scant 21 km from Lituya Bay. Intense shaking lasted for 1 to 4 min according to two eyewitnesses who were anchored at the entrance of the bay. According to Miller (1960), between 1 and 2.5 min after the earthquake, a large mass of rock slid from the northeast wall of Gilbert Inlet. It is probable that this entire mass of rocks, ice, and soil plunged into Gilbert Inlet as a unit. The result was the sudden displacement of a large volume of water as the slide was plunged into Gilbert Inlet, causing the largest tsunami ever evidenced. The upper limit of destruction of forest and vegetation by water (known as trimline) extended to a maximum of 524 m a.m.s.l. (above mean sea level) on the spur southwest of Gilbert Inlet (Fig. 1). Maximum inundation distance reached 1400 m on flat ground at Fish Lake on the north side of the bay, closer to its entrance.

In order to understand the evolution of the giant Lituya Bay wave, a rough model at a 1 : 1000 scale was constructed at the University of California, Berkeley (Robert L. Wiegel in Miller, 1960). If the slide occurred rapidly as a unit, they concluded that a sheet of water washed up the slope opposite the landslide to an elevation of at least 3 times the water depth. At the same time, a large wave, several hundred feet high, moved in a southward direction, causing a peak rise to occur in the vicinity of Mudslide Creek, on the south shore of Gilbert Inlet. According to Miller (1960), this peak reached 204 m (580 ft) (see Fig. 1).

The landslide was triggered by fault movement and intense earthquake vibrations (Fritz et al., 2009) and it is highly probable that the entire mass of rock plunged into Gilbert Inlet as a unit, as previously stated. Nevertheless, there is no consensus about the typology of the slide mass movement. Miller (1960) provides discussion setting this event near the borderline between a landslide and a rockfall following the classifications of Sharpe (1938) and Varnes (1958), while Pararas-Carayannis (1999) classified the mass movement as subaerial rockfall, making the distinction from gradual processes of ordinary landslides. Nevertheless, as will be shown in the next sections, in Fritz et al. (2001, 2009), the authors proposed, and in fact show a landslide typology that, based on the generalized Froude similarity, makes it possible to reproduce this event using a two-dimensional scaled physical model of the Gilbert Inlet. A pneumatic landslide generator was used to generate a high-speed granular slide with density and volume based on Miller (1960) impacting the water surface at a mean velocity of 110 m s⁻¹. The experimental run-up matches the trimline of forest destruction on the spur ridge in Gilbert Inlet.

Figure 1 Topographic map of Lituya Bay (US Geologic Survey, 1961) showing the settings and trimline of the 1958 mega-tsunami (based on data from Miller, 1960). Units in feet. Key locations such as Gilbert and Crillon inlets, Cenotaph Island, Fish Lake, or The Paps are shown.

Based on the experimental work of Fritz et al. (2009), in the present numerical study we will follow the same approach: an initial slide speed (analogous to the impulse of the pneumatic landslide generator in the lab experiment) will be imposed in order to obtain the 110 m s⁻¹ slide impact velocity that Fritz et al. (2009) measure in their experiment. In the same way that the laboratory experiment reproduced the observed run-up, it is the way in which the numerical experiments performed here have been initialized.

3 Area of study

Lituya Bay (Fig. 1), located within Glacier Bay National Park on the northeast shore of the Gulf of Alaska, is a T-shaped tidal inlet, nearly 12 km long and with a width ranging from 1.2 to 3.3 km except at the entrance, where the width is approximately 300 m. The northeastward stem of the bay cuts the coastal lowlands and the foothills flanking the Fairweather Range of the St. Elias Mountains. In the vicinity of the head of the bay, the walls are steep, fjord-like, and rise to elevations ranging from 670 to 1040 m in the foothills immediately to the north and south, and more than 1800 m in the Fairweather Range. In 1958 the maximum depth of the bay was 220 m and the sill depth, at the entrance of the bay, was only 10 m. At the head of the bay, the two arms of the T, the Gilbert (northern arm) and Crillon (southern arm) inlets, form part of a great trench that extends tens of kilometers to the northwest and southeast on the Fairweather Fault. Cenotaph Island divides the central part of the bay into two channels of 640 and 1290 m, respectively.

Figure 2 Gilbert Inlet shoreline and glacier, bathymetries, and shorelines before and after 1958. (a) Modified location of shoreline and glacier front before 1958, from Miller (1960). (b) Shoreline and glacier front of Lituya Glacier after the 1958 tsunami. (c) Reconstructed Gilbert Inlet bathymetry (based on 1926 data). (d) Gilbert Inlet bathymetry (1959).

3.1 Coastal morphology

The shores around the main part of the bay are composed mainly of rocky beaches that rise steeply away from the shoreline. There are two adjoining land masses that rise away from the beach, ranging in elevations from less than 30 m at a horizontal distance of 2 km, around Fish Lake, to 170 m at a horizontal distance of 370 m at The Paps (see Fig. 1). Prior to the 1958 tsunami, low deltas of gravel had built out into Gilbert Inlet along the southwest and northeast margins of the Lituya Glacier front.

According to Miller (1960), and as evidenced in several graphical documents, after the tsunami, the delta on the northeast side of Gilbert Inlet completely disappeared, and that on the southwest side of the bay was noted to be significantly smaller. To recreate the scenario previous to the event for use in the numerical simulation presented here, shorelines, deltas, and the glacial front inside Gilbert Inlet before 1958 were all taken from Miller's reconstruction (see Fig. 2).

4 Bathymetric data

According to Miller (1960), examination of Lituya Bay bathymetry was the first step in determining whether the volume of water was sufficient to account for a 524 m wave. Bathymetric surveys made in 1926 and 1940 (US Coast and Geodesic Survey, 1926) show that the head of Lituya Bay is a pronounced U-shaped trench with steep walls and a broad, flat floor that slopes gently downward from the head of the bay to a maximum depth of 220 m just south of Cenotaph Island. At this maximum depth, the slope then rises toward the outer part of the bay. At the entrance of the bay, the minimum depth is on the order of 10 m during mean lower low water. The outer portion of Lituya Bay is enclosed by a long spit, the La Chaussee Spit, with only a very narrow entrance of about 220–245 m that is kept open by tidal currents. The tide in the bay is predominantly diurnal, with a mean range of 2 m and a maximum range of about 4.5 m (US Coast and Geodesic Survey, 1959). The U shape of the bay coupled with the flatness of its floor indicate that extensive sedimentation has taken place. The thickness of the sediments at the Gilbert and Crillon inlets is not known but believed to be substantial due to terminal moraine deposition during different brief glacial and interglacial episodes.

The bathymetric data used for the modeling work presented here were obtained from the US National Ocean Service: Hydrographic Surveys with Digital Sounding. Data from survey ID H08492, 1959, were used as reference bathymetry since this survey is the nearest in time to the data of the tsunami and there were enough data collected to provide a good representation of the entire Lituya Bay seafloor. Data from survey ID H04608, 1926, were used to reconstruct Gilbert Inlet bathymetry as these were the closest pre-event data available. Unfortunately, data from this survey are not sufficient in resolution to provide an acceptable bathymetric grid for our study of the entire bay. Nevertheless, the survey provides both enough data and detailed information of pre-tsunami bathymetry in the Gilbert Inlet. In Fig. 2c and d the 1926 reconstructed bathymetry and the original 1959 bathymetry are shown, respectively. The mass volume difference between these two bathymetries in the Gilbert Inlet area is about 31× 10⁶ m³, that is, very close to the slide volume estimated by Miller (1960).

5 Tsunami source

The dimensions of the landslide on the northeast wall of Gilbert Inlet were determined with reasonable accuracy by Miller (1960), but the thickness of the slide mass normal to the slope could be estimated only roughly from available data and photographs. The main mass of the slide was a prism of rock that was roughly triangular in cross section, with dimensions from 730 to 915 m along the slope, a maximum thickness of about 92 m normal to the slope, and a center of gravity at about 610 m in elevation. From these dimensions the volume of the slide estimated by Miller (1960) was of 30.6 × 10⁶ m³.

To locate and reconstruct the volume of the slide mass, the following procedure was implemented: first, based on aerial photos and data provided by Miller (1960), the perimeter of the slide was determined. Then, an approximate centroid for the formerly defined surface was considered drawing two lines, one horizontal and another one vertically projected on the surface. The surface centroid was located at 610 m in height, defining the upper bound for the mass slide. The volume of the reconstructed slide was 30.625×10⁶ m³, which matches accurately with Miller's estimation.

The three criteria we tried to fulfill in order to reconstruct the slide geometry were (1) to place it in its exact location (projected area), (2) keep an approximate location for the centroid, also in height, and (3) to recover an accurate volume for the numerical slide.

5.1 Landslide setup

In order to reproduce the main features of the slide impact, Hermann M. Fritz and collaborators (Fritz et al., 2009) designed a pneumatic landslide generator. They intend to model the transition from rigid to granular slide motion. Thus, at the beginning, the granular material is impulsed until the landslide achieves 110 m s⁻¹, which is the approximated impact velocity between the slide and the water surface estimated by Fritz et al., assuming free-fall equations for the centroid of the slide. From this instant, the slide is supposed to behave as a granular medium.

In this work, we have followed the same idea: assuming that 110 m s⁻¹ is a good approximation for the impact velocity of the slide, an initial velocity for the granular layer has been estimated so that the computed impact velocity is approximately 110 m s⁻¹. This is a critical point for the performance of the simulation in order to reproduce the dynamic of the impact, the generation of the tsunami, the propagation, and the run-up on different areas of the domain. Thus, for example, if model simulation is initialized from a solid slide at rest, this results in an impact velocity of the slide with the water of only, approximately, 67 m s⁻¹.

From a detailed analysis of the bathymetry surveys available for this study, an unexpected shifted location of the slide deposit on the floor of the Gilbert Inlet (see Fig. 2c), with a larger deposit concentration to the south part of the inlet, was observed. The observation of this fact made us consider a slide initial velocity vector, v_s, slightly shifted to the south, with a modulus close to 80 m s⁻¹ (see Fig. 3). Thus, with this initial condition, the model reproduces both the run-up on the spur southwest of Gilbert Inlet and a giant wave traveling into the bay with enough energy to accurately reproduce the effects of the wave along the Lituya Bay.

6 Model description

Coulomb-type models for granular-driven flows have been intensively investigated in the last decade, following the pioneering work of Savage and Hutter (Savage and Hutter, 1989), who derived a shallow-water-type model including a Coulomb friction term to take into account the interaction of the avalanche with the bottom topography. This model has been extended and generalized in several works (Bouchut et al., 2003, 2008, 2016; Bouchut and Westdickenberg, 2004; Pelanti et al, 2008). In this framework, the EDANYA group has implemented a finite-volume numerical model, the Landslide-HySEA model (Castro et al., 2005, 2006, 2008, 2012; Gallardo et al., 2007; Macías et al., 2012, 2015; de la Asunción et al., 2013), for the simulation of submarine landslides based on the two-layer Savage–Hutter model introduced in Fernández-Nieto et al. (2008) that takes into account the movement of the fluid inside which the avalanche develops. This model is useful for the generation and evolution of tsunamis triggered by both submarine and aerial landslides.

HySEA models have been fully validated for tsunami modeling using all of the benchmark problems posted on the NOAA NTHMP web site for propagation and inundation (Macías et al., 2017b), for tsunami currents (Macías et al., 2019a, b) and for landslides (Macías et al., 2017a). In Sánchez-Linares (2011), the Landslide-HySEA model is used to reproduce the Fritz et al. (2009) laboratory experiment, producing good results.

Figure 3 Initial conditions for the slide: slide location and initial velocity vector direction.

This section describes the system of partial differential equations modeling landslide-generated tsunamis based on layered average models. The 2-D Landslide-HySEA model is a two-dimensional version of the model proposed in Fernández-Nieto et al. (2008) for 1-D problems, for which local coordinates are not considered.

6.1 Simplified two-layer Savage–Hutter-type model

Let us consider a layered medium composed by a layer of inviscid fluid with constant homogeneous density ρ₁ (water) and a layer of granular material with density ρ_s and porosity ψ₀. We assume that both layers are immiscible and the mean density of the granular material layer is given by ${\mathit{\rho }}_{\mathrm{2}}=(\mathrm{1}-{\mathit{\psi }}_{\mathrm{0}}){\mathit{\rho }}_{\mathrm{s}}+{\mathit{\psi }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{1}}$. The system of partial differential equations (PDEs) describing the coupled two-layer system is written as

$$\begin{array}{}\text{(1)}& {\displaystyle }\left\{\begin{array}{l}{\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial t}}+{\displaystyle \frac{\partial q_{\mathrm{1},x}}{\partial x}}+{\displaystyle \frac{\partial q_{\mathrm{1},y}}{\partial y}}=\mathrm{0}\\ \\ {\displaystyle \frac{\partial q_{\mathrm{1},x}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{q_{\mathrm{1},x}^{\mathrm{2}}}{h_{\mathrm{1}}}}+{\displaystyle \frac{g}{\mathrm{2}}}{h}_{\mathrm{1}}^{\mathrm{2}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{q_{\mathrm{1},x}{q}_{\mathrm{1},y}}{h_{\mathrm{1}}}}\right)=\\ -g{h}_{\mathrm{1}}{\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial x}}+g{h}_{\mathrm{1}}{\displaystyle \frac{\partial H}{\partial x}}+S_{{\mathrm{f}}_{\mathrm{1}}}\left(\mathit{W}\right)\\ \\ {\displaystyle \frac{\partial q_{\mathrm{1},y}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{q_{\mathrm{1},x}{q}_{\mathrm{1},y}}{h_{\mathrm{1}}}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{q_{\mathrm{1},y}^{\mathrm{2}}}{h_{\mathrm{1}}}}+{\displaystyle \frac{g}{\mathrm{2}}}{h}_{\mathrm{1}}^{\mathrm{2}}\right)=\\ -g{h}_{\mathrm{1}}{\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial y}}+g{h}_{\mathrm{1}}{\displaystyle \frac{\partial H}{\partial y}}+S_{{\mathrm{f}}_{\mathrm{2}}}\left(\mathit{W}\right)\\ \\ {\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial t}}+{\displaystyle \frac{\partial q_{\mathrm{2},x}}{\partial x}}+{\displaystyle \frac{\partial q_{\mathrm{2},y}}{\partial y}}=\mathrm{0}\\ \\ {\displaystyle \frac{\partial q_{\mathrm{2},x}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{q_{\mathrm{2},x}^{\mathrm{2}}}{h_{\mathrm{2}}}}+{\displaystyle \frac{g}{\mathrm{2}}}{h}_{\mathrm{2}}^{\mathrm{2}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{q_{\mathrm{2},x}{q}_{\mathrm{2},y}}{h_{\mathrm{2}}}}\right)=\\ -gr{h}_{\mathrm{2}}{\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial x}}+g{h}_{\mathrm{2}}{\displaystyle \frac{\partial H}{\partial x}}+S_{{\mathrm{f}}_{\mathrm{3}}}\left(\mathit{W}\right)+{\mathit{\tau }}_x\\ \\ {\displaystyle \frac{\partial q_{\mathrm{2},y}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{q_{\mathrm{2},x}{q}_{\mathrm{2},y}}{h_{\mathrm{2}}}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{q_{\mathrm{2},y}^{\mathrm{2}}}{h_{\mathrm{2}}}}+{\displaystyle \frac{g}{\mathrm{2}}}{h}_{\mathrm{2}}^{\mathrm{2}}\right)=\\ -gr{h}_{\mathrm{2}}{\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial y}}+g{h}_{\mathrm{2}}{\displaystyle \frac{\partial H}{\partial y}}+S_{{\mathrm{f}}_{\mathrm{4}}}\left(\mathit{W}\right)+{\mathit{\tau }}_y.\end{array}\right.\end{array}$$

Figure 4 Details about the reconstructed slide. The isolines into the yellow area have been modified to reconstruct the location of the slide previous to the event. (a) Reconstructed topography of the east spur of Gilbert Inlet (before). (b) Topography of the east spur of Gilbert Inlet after the tsunami.

In these equations, the subscript 1 refers to the fluid upper layer, and subscript 2 to the lower layer composed of the fluidized material. The ith layer thickness at point $(x,y)\in D\subset {\mathbb{R}}^{\mathrm{2}}$ at time t, where D is the horizontal projection of domain occupied by the fluid, is denoted by $h_i(x,y,t)$. H(x,y) indicates the depth of non-erodible bottom measured from a fixed reference level at point (x,y), and ${\mathit{q}}_i(x,y,t)=\left(q_{i,x}(x,y,t),q_{i,y}(x,y,t)\right)$ is the flow of the ith layer at point (x,y) at time t, which are related to the mean velocity of each layer (${\mathit{u}}_i(x,y,t)$) by ${\mathit{q}}_i(x,y,t)=h_i(x,y,t){\mathit{u}}_i(x,y,t)$, $i=\mathrm{1},\mathrm{2}$. The value $r={\mathit{\rho }}_{\mathrm{1}}/{\mathit{\rho }}_{\mathrm{2}}$ denotes the ratio between the constant densities of the two layers (ρ₁<ρ₂). Note that H(x,y) does not depend on t, that is, the non-erodible bottom topography does not change through the simulation although the bottom may change due to second layer movement. Figure 5 graphically shows the relationship among h₁, h₂, and H. Usually, $h_{\mathrm{1}}+h_{\mathrm{2}}=H$ is at rest or it represents the mean sea level.

The terms $S_{{\mathrm{f}}_i}\left(\mathit{W}\right)$, $i=\mathrm{1},\mathrm{\dots },\mathrm{4}$ (where W refers to model variables, i.e., $\mathit{W}=(h_{\mathrm{1}};q_{\mathrm{1},x};q_{\mathrm{1},y};h_{\mathrm{2}};q_{\mathrm{2},x};q_{\mathrm{2},y})$) model the different effects of dynamical friction, while $\mathit{\tau }=({\mathit{\tau }}_x,{\mathit{\tau }}_y)$ is the Coulomb static friction law. $S_{{\mathrm{f}}_i}\left(\mathit{W}\right)$, $i=\mathrm{1},\mathrm{\dots },\mathrm{4}$ are given by

$$\begin{array}{}\text{(2)}& {\displaystyle }\left\{\begin{array}{l}{S}_{{\mathrm{f}}_{\mathrm{1}}}\left(\mathit{W}\right)=S_{c_x}\left(\mathit{W}\right)+S_{a_x}\left(\mathit{W}\right);\\ S_{{\mathrm{f}}_{\mathrm{2}}}\left(\mathit{W}\right)=S_{c_y}\left(\mathit{W}\right)+S_{a_y}\left(\mathit{W}\right);\\ \\ S_{{\mathrm{f}}_{\mathrm{3}}}\left(\mathit{W}\right)=-r{S}_{c_x}\left(\mathit{W}\right)+S_{{\mathrm{b}}_x}\left(\mathit{W}\right);\\ S_{{\mathrm{f}}_{\mathrm{4}}}\left(\mathit{W}\right)=-r{S}_{c_y}\left(\mathit{W}\right)+S_{{\mathrm{b}}_y}\left(\mathit{W}\right).\end{array}\right.\end{array}$$

The term ${\mathit{S}}_c\left(\mathit{W}\right)=\left(S_{c_x}\right(\mathit{W}),S_{c_y}(\mathit{W}\left)\right)$ parameterizes the friction between the two layers, water and granular slide, as a quadratic function of the shear velocity and is defined (as a particular case of the parameterization proposed in Pitman and Le, 2005, or Pelanti et al, 2008) as

$$\begin{array}{}\text{(3)}& {\displaystyle }\left\{\begin{array}{l}{S}_{c_x}\left(\mathit{W}\right)=m_{\mathrm{f}}{\displaystyle \frac{h_{\mathrm{1}}{h}_{\mathrm{2}}}{h_{\mathrm{2}}+r{h}_{\mathrm{1}}}}(u_{\mathrm{2},x}-u_{\mathrm{1},x})\Vert {\mathit{u}}_{\mathrm{2}}-{\mathit{u}}_{\mathrm{1}}\Vert \\ \\ S_{c_y}\left(\mathit{W}\right)=m_{\mathrm{f}}{\displaystyle \frac{h_{\mathrm{1}}{h}_{\mathrm{2}}}{h_{\mathrm{2}}+r{h}_{\mathrm{1}}}}(u_{\mathrm{2},y}-u_{\mathrm{1},y})\Vert {\mathit{u}}_{\mathrm{2}}-{\mathit{u}}_{\mathrm{1}}\Vert ,\end{array}\right.\end{array}$$

where m_f is a positive constant, and ${\mathit{S}}_a\left(\mathit{W}\right)=\left(S_{a_x}\right(\mathit{W}),S_{a_y}(\mathit{W}\left)\right)$ parameterizes the friction between the fluid and the non-erodible bottom and is given by a Manning law (Dyakonova and Khoperskov, 2018):

$$\begin{array}{}\text{(4)}& {\displaystyle }\left\{\begin{array}{l}{S}_{a_x}\left(\mathit{W}\right)=-g{h}_{\mathrm{1}}{\displaystyle \frac{n_{\mathrm{1}}^{\mathrm{2}}}{h_{\mathrm{1}}^{\mathrm{4}/\mathrm{3}}}}{u}_{\mathrm{1},x}\Vert {\mathit{u}}_{\mathrm{1}}\Vert \\ S_{a_y}\left(\mathit{W}\right)=-g{h}_{\mathrm{1}}{\displaystyle \frac{n_{\mathrm{1}}^{\mathrm{2}}}{h_{\mathrm{1}}^{\mathrm{4}/\mathrm{3}}}}{u}_{\mathrm{1},y}\Vert {\mathit{u}}_{\mathrm{1}}\Vert ,\end{array}\right.\end{array}$$

where n₁>0 is the Manning coefficient.

${\mathit{S}}_{\mathrm{b}}\left(\mathit{W}\right)=\left(S_{{\mathrm{b}}_x}\right(\mathit{W}),S_{{\mathrm{b}}_y}(\mathit{W}\left)\right)$ parameterizes the friction between the granular and the non-erodible bottom and as in the previous case is given by a Manning law:

$$\begin{array}{}\text{(5)}& {\displaystyle }\left\{\begin{array}{l}{S}_{{\mathrm{b}}_x}\left(\mathit{W}\right)=-g{h}_{\mathrm{2}}{\displaystyle \frac{n_{\mathrm{2}}^{\mathrm{2}}}{h_{\mathrm{2}}^{\mathrm{4}/\mathrm{3}}}}{u}_{\mathrm{2},x}\Vert {\mathit{u}}_{\mathrm{2}}\Vert \\ S_{{\mathrm{b}}_y}\left(\mathit{W}\right)=-g{h}_{\mathrm{2}}{\displaystyle \frac{n_{\mathrm{2}}^{\mathrm{2}}}{h_{\mathrm{2}}^{\mathrm{4}/\mathrm{3}}}}{u}_{\mathrm{2},y}\Vert {\mathit{u}}_{\mathrm{2}}\Vert ,\end{array}\right.\end{array}$$

where n₂>0 is the corresponding Manning coefficient.

Note that S_a(W) is only defined where $h_{\mathrm{2}}(x,y,t)=\mathrm{0}$. In this case, m_f=0 and n₂=0. Similarly, if $h_{\mathrm{1}}(x,y,t)=\mathrm{0}$ then m_f=0 and n₁=0.

Finally, $\mathit{\tau }=({\mathit{\tau }}_x,{\mathit{\tau }}_y)$ is defined as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{l}\text{If }\Vert \mathit{\tau }\Vert \ge {\mathit{\sigma }}^c\Rightarrow \left\{\begin{array}{l}{\mathit{\tau }}_x=-g(\mathrm{1}-r)h_{\mathrm{2}}{\displaystyle \frac{q_{\mathrm{2},x}}{\Vert {\mathit{q}}_{\mathrm{2}}\Vert }}\tan \left(\mathit{\alpha }\right)\\ {\mathit{\tau }}_y=-g(\mathrm{1}-r)h_{\mathrm{2}}{\displaystyle \frac{q_{\mathrm{2},y}}{\Vert {\mathit{q}}_{\mathrm{2}}\Vert }}\tan \left(\mathit{\alpha }\right)\end{array}\right.\\ \text{If }\Vert \mathit{\tau }\Vert <{\mathit{\sigma }}^c\Rightarrow q_{\mathrm{2},x}=\mathrm{0},q_{\mathrm{2},y}=\mathrm{0}\end{array},\end{array}$$

where ${\mathit{\sigma }}^c=g(\mathrm{1}-r)h_{\mathrm{2}}\tan \left(\mathit{\alpha }\right)$, where α is the Coulomb friction angle (Savage and Hutter, 1989; Gray et al., 1999).

System Eq. (1) can be written as a system of conservation laws with source terms and nonconservative products (Fernández-Nieto et al., 2008). In the next section, the finite-volume scheme used to discretize system Eq. (1) is described. As friction terms are semi-implicitly discretized, we first consider that S_F(W)=0, with S_F being the vector containing all friction term contributions, i.e., ${\mathit{S}}_{\mathrm{F}}\left(\mathit{W}\right)=\left(S_{{\mathrm{f}}_{\mathrm{1}}}\right(\mathit{W});S_{{\mathrm{f}}_{\mathrm{2}}}(\mathit{W});S_{{\mathrm{f}}_{\mathrm{3}}}(\mathit{W});S_{{\mathrm{f}}_{\mathrm{4}}}(\mathit{W}\left)\right)$. Then, the way those terms are discretized is briefly described.

7 Numerical scheme

To discretize system Eq. (1), the domain D is divided into L cells or finite volumes V_i⊂ℝ², $i=\mathrm{1},\mathrm{\dots },L$, which are assumed to be closed polygons. We assume here that the cells are rectangles with edges parallel to Cartesian axes. Given a finite volume V_i, N_i∈ℝ² is the center of V_i, ℵ_i is the set of indices j such that V_j is a neighbor of V_i, Γ_ij is the common edge of two neighboring volumes V_i and V_j, and $\left|{\mathrm{\Gamma }}_{ij}\right|$ is its length; ${\mathit{\eta }}_{ij}=({\mathit{\eta }}_{ij,x},{\mathit{\eta }}_{ij,y})$ is the unit normal vector to the edge Γ_ij and pointing towards V_j.

We denote by ${\mathit{W}}_i^n$ an approximation of the solution average over the cell V_i at time tⁿ:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{W}}_i^n\cong {\displaystyle \frac{\mathrm{1}}{\left|V_i\right|}}\underset{V_i}{\int }\mathit{W}(x,y,t^n)\mathrm{d}x\mathrm{d}y,\end{array}$$

where $\left|V_i\right|$ is the area of cell V_i and $t^n=t^{n-\mathrm{1}}+\mathrm{\Delta }t$, where Δt is the time step.

Let us suppose that ${\mathit{W}}_i^n$ is known. Thus, to advance in time, a family of one-dimensional Riemann problems projected in the normal direction to each edge of the mesh Γ_ij is considered. Those Riemann problems are approximated by means of an intermediate field capturing parabola (IFCP) numerical scheme (see Fernández-Nieto et al., 2011). Finally, ${\mathit{W}}_i^{n+\mathrm{1}}$ is computed by averaging these approximate solutions. The resulting numerical scheme is written as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{W}}_i^{n+\mathrm{1}}={\mathit{W}}_i^n-{\displaystyle \frac{\mathrm{\Delta }t}{\left|V_i\right|}}\sum _{j\in {\mathrm{\aleph }}_i}\left|{\mathrm{\Gamma }}_{ij}\right|{\mathcal{F}}_{ij}^{-}({\mathit{W}}_i^n,{\mathit{W}}_j^n,H_i,H_j).\end{array}$$

To check the precise definition for the numerical fluxes, ${\mathcal{F}}_{ij}^{-}({\mathit{W}}_i^n,{\mathit{W}}_j^n,H_i,H_j)$, see Sánchez-Linares et al. (2015) or de la Asunción et al. (2016).

7.1 Wet–dry fronts

The numerical scheme described above, when applied to wet–dry situations, may produce incorrect results: the gradient of the bottom topography generates spurious pressure forces and the fluid can artificially climb up slopes. In Castro et al. (2005), to avoid this problem, a modification of the numerical scheme is proposed. Here the same strategy is to correct the proposed numerical scheme to suitably deal with wet–dry fronts. With this strategy, spurious wave reflections in the coast are avoided and a more realistic simulation of the flooded areas is obtained. Moreover, transitions between sub- and supercritical flows, which appear continuously in simulations such as the one presented here, which further complicate matters, are also suitably treated.

The implementation of the wet–dry front treatment in the numerical scheme results in not having to impose boundary conditions at the coasts. The coastline becomes a moving boundary, computed by the numerical scheme. Depending on the impact wave characteristics or the water backflow movement, the computational cells are filled with water or they run dry. Consequently, no specific stabilization model technique is required.

7.2 Friction term discretization

In this section the numerical scheme when the friction terms S_F(W) are discretized is presented. First, the terms $S_{{\mathrm{f}}_{\mathrm{1}}}\left(\mathit{W}\right)$, $S_{{\mathrm{f}}_{\mathrm{2}}}\left(\mathit{W}\right)$, $S_{{\mathrm{f}}_{\mathrm{3}}}\left(\mathit{W}\right)$, and $S_{{\mathrm{f}}_{\mathrm{4}}}\left(\mathit{W}\right)$ are discretized semi-implicitly; next the Coulomb friction term τ will be discretized following Fernández-Nieto et al. (2008). The resulting numerical scheme is a three-step method, in which in the first step the IFCP scheme is used and then in the other two steps the dynamical and static friction terms will be discretized. The three-step method will be denoted as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{W}}_i^n\to {\mathit{W}}_i^{n+\mathrm{1}/\mathrm{3}}\to {\mathit{W}}_i^{n+\mathrm{2}/\mathrm{3}}\to {\mathit{W}}_i^{n+\mathrm{1}}.\end{array}$$

Figure 5 Sketch of the two-layer model. Relation among h₁, h₂, and H.

The resulting scheme is exactly well balanced (Sánchez-Linares et al., 2015) for the stationary water at rest solution (${\mathit{q}}_{\mathrm{1}}={\mathit{q}}_{\mathrm{2}}=\mathbf{0}$ and μ₁ and μ₂ constant). Moreover, the scheme accurately solves the stationary solutions corresponding to ${\mathit{q}}_{\mathrm{1}}={\mathit{q}}_{\mathrm{2}}=\mathbf{0}$, μ₁ constant and ${\partial }_x{\mathit{\mu }}_{\mathrm{2}}<\tan \left(\mathit{\alpha }\right)$ and ${\partial }_y{\mathit{\mu }}_{\mathrm{2}}<\tan \left(\mathit{\alpha }\right)$, that is, a stationary water at rest solution for which the Coulomb friction term balances the pressure term in the granular material.

7.3 Numerical resolution and GPU implementation

Landslide initial conditions have been described in Sect. 5.1. Initially water is at rest, and according to Miller (1960), an initial level of 1.52 m has been set.

The computational grid considered for the numerical simulation is a 4 m × 7.5 m rectangular mesh composed of 3650×1271, i.e., 4 639 150, cells.

The numerical scheme described in this section has been implemented in the CUDA programming language in order to be able to run the model in GPUs. A highly efficient GPU-based implementation of the numerical scheme allows us to compute highly accurate simulations in very reasonable computational time.

The numerical simulation presented in this work covers a wall clock time period of 10 min. A total of 14 516 time iterations were necessary to evolve from initial conditions to final state, 10 min later. This required a computational time of 1528.83 s (approx 25.5 min), which means 44 million computational cells were processed per second in a nVidia GTX 480 graphic card.

8 Parameter sensitivity analysis

To overcome the uncertainty inherent in the choice of model parameters and in order to produce a numerical simulation as close as possible to the real event, a sensitivity analysis has been performed. To do so, the three key parameters, (1) Coulomb friction angle, α, (2) the ratio of densities between the water and the mean density of the slide, r, and (3) the friction between layers, m_f, have been retained as varying parameters for this sensitivity analysis. The values for these three parameters have been moved over the following ranges of reasonable values:

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathit{\alpha }\in [\mathrm{8}{}^{\circ },\mathrm{16}{}^{\circ }]r\in [\mathrm{0.3},\mathrm{0.6}]m_{\mathrm{f}}\in [\mathrm{0.001},\mathrm{0.1}].\end{array}$$

The other two parameters required in model parameterizations, the Manning coefficients n₁ and n₂, were set constant with standard values of n₁=0.02 and n₂=0.05 (Arcement and Schneider, 1989; Phillips and Tadayon, 2006).

Figure 6 Location and spatial extension of the four regions considered to compute maximum run-up for the sensitivity analysis.

Previous to searching the optimal set of parameters that will determine the optimal solution, a sensitivity analysis was carried out. In order to assess model sensitivity to parameters, four regions have been selected (see Fig. 6) to measure the maximum run-up at the coastal strip within these regions. These values will serve as a proxy for model sensitivity and will be compared with the observed values provided by Miller (1960). To do this, a “macroscopic” (in the sense that a larger sampling step is used) set of parameters is considered. Then, the optimal solution will be searched using a local refinement of the former sampling in order to obtain more precise values for the final optimal parameters. For the macroscopic set of parameters the following values were used: for r, 0.3, 0.4, 0.5, and 0.6; for α, 8, 10, 12, 14, and 16^∘; and for m_f, 0.001, 0.005, 0.01, 0.025, 0.05, 0.075, and 0.1. Summing up a total of 140 numerical experiments. The numerical results have been summarized in 3-by-4 graphs presented here as three figures (Figs. 7, 8, and 9).

Figure 7 For a fixed value of the Coulomb friction angle, α (in (a) $\mathit{\alpha }=\mathrm{8}{}^{\circ }$; (b) $\mathit{\alpha }=\mathrm{10}{}^{\circ }$; (c) $\mathit{\alpha }=\mathrm{14}{}^{\circ }$; and (d) $\mathit{\alpha }=\mathrm{16}{}^{\circ }$), the maximum computed run-up in the four regions defined in Fig. 6 for different values of the friction coefficient, m_f (in different colors; see legend), is presented. The vertical bars show the run-up variation for values of r ranging from 0.3 to 0.6. The horizontal lines (in black) show the reference values for each region: (A) 524 m, (B) 165 (solid) and 122 m (dashed), (C) higher than 60 m, and (D) 50 m. In addition, horizontal lines in red show the values for the optimal solution that will be described below.

Each panel in Fig. 7 summarizes the variability in the maximum computed run-up in 35 (out of the 140) of the numerical experiments composing the macroscopic sampling, in each of the four regions defined in Fig. 6 for a value of the Coulomb friction angle, α, fixed: $\mathit{\alpha }=\mathrm{8}{}^{\circ }$ in (A), $\mathit{\alpha }=\mathrm{10}{}^{\circ }$ in (B), $\mathit{\alpha }=\mathrm{14}{}^{\circ }$ in (C), and $\mathit{\alpha }=\mathrm{16}{}^{\circ }$ in (D). In different colors, as indicated by the legend, the maximum run-up for different values of the friction coefficient, m_f, is presented. Finally, the vertical bars show the run-up variation for values of r ranging from 0.3 to 0.6 for the same friction coefficient. The reference values, as extracted from Miller (1960), are also presented (as horizontal black lines) in all figures, and these are as follows. For region (A) the maximum run-up is 524 m. For region (B) the two maximum values (165 and 122 m) indicated by Miller (1960) in this area have been marked (it was noted that in this area the model underestimates Miller's observed values). For region (C), for water to cross the Cenotaph Island, a run-up larger than 60 m is required. For region (D) 50 m (see Fig. 10d) is set as reference. The horizontal red lines show the values for the optimal solution that will be determined below.

From this first set of figures, we can conclude that m_f=0.001 it is not an admissible value for the interface friction and that m_f=0.05, although already in the range of observed values for most cases, still seems not to be suitable. This is the reason why in the next set of figures (in order to restrict ourselves to four panels), the figure for m_f=0.001 has not been included. The figures for m_f=0.05 and m_f=0.075 are not presented since the results are well represented by the neighbor values m_f=0.025 and m_f=0.1.

In Fig. 8, each panel shows the values of the maximum run-up computed in each region for a fixed value of the friction coefficient, m_f, in particular m_f=0.005 for (A), m_f=0.01 in (B), m_f=0.025 in (C), and m_f=0.1 in (D). The colored lines gather the run-ups for different values of the Coulomb friction angle, α. Finally, the vertical bars show the run-up variation for values of r ranging from 0.3 to 0.6. Again, it can be seen that the smaller values for m_f produce numerical values farther from observations, and values of the friction between layers between 0.025 and 0.1 seem to produce numerical results closer to observations.

Figure 8 For a fixed value of the friction coefficient, m_f (in (a) m_f=0.005; (b) m_f=0.01; (c) m_f=0.025; and (d) m_f=0.1), the maximum computed run-up in the four regions defined in Fig. 6 for different values of the Coulomb friction angle, α (in different colors; see legend), is presented. The vertical bars show the run-up variation for values of r ranging from 0.3 to 0.6.

In Fig. 9, each panel shows the values of the maximum run-up computed in each region for a fixed value of the ratio of densities, r, in particular r=0.3 in (A), r=0.4 in (B), r=0.5 in (C), and r=0.6 in (D). The colored lines gather the run-ups for different values of the friction coefficient, m_f. Finally, the vertical bars show the run-up variation for values of α ranging from 8 to 16^∘. This set of figures clearly points out that the value for r must be between 0.4 and 0.5, and as in previous figures, it can be observed that the smaller values for the friction parameter must be ruled out.

Figure 9 For a fixed value of the ratio of densities, r (in (a) r=0.3; (b) r=0.4; (c) r=0.5; and (d) r=0.6), the maximum computed run-up in the four regions defined in Fig. 6 for different values of the friction coefficient, m_f (in different colors; see legend), is presented. The vertical bars show the run-up variation for values of α ranging from 8 to 16^∘.

Figure 10 Images showing the four criteria used to select the optimal numerical solution and comparison of simulated optimal inundated area with the observed trimline. (a) Gilbert and Crillon inlets. (b) Inner part of the Lituya Bay from Cenotaph Island. (c) Fish Lake and Cenotaph Island area. (d) Entrance of the bay up to Cenotaph Island.

Now, in order to find the optimal parameters, according to four criteria that will be defined below, a more detailed “microscopic” grid of parameters  is used for performing additional numerical experiments. Taking into account the results presented in Figs. 7, 8, and 9, the ratio of densities, r, will now be varied between 0.4 and 0.5, the Coulomb friction angle, α, between 12 and 14^∘, and the friction between layers, m_f, between 0.01 and 0.1. More specifically, for this microscopic set of parameters, the following values were used: for r, $\mathrm{0.4},\mathrm{0.42},\mathrm{0.44},\mathrm{0.46},\mathrm{0.48}$, and 0.5; for α, 12, 13, and 14^∘; and for m_f, 0.01, 0.02, 0.04, 0.06, 0.08, and 0.1. Summing up a total of 108 additional numerical experiments.

Figure 11 Giant wave generation and initial evolution in Gilbert Inlet. (a) t=0 s; (b) t=8 s. The giant wave reaches 272.4 m in amplitude; (c) t=10 s. Maximum wave amplitude reaches 251.1 m. (d) t=20 s. Maximum wave amplitude reaches 161.5 m; (e) t=30 s. The giant wave hits the spur southwest of Gilbert Inlet; (f) t=39 s. Maximum run-up: 523.9 m.

Then, the four criteria considered in order to select the optimal parameters were as follows.

C1

The run-up on the spur southwest of Gilbert Inlet had to be the closest to the optimal 524 m.

C2

The wave moving southwards to the main stem of Lituya bay had to cause a peak close to 208 m in the vicinity of Mudslide Creek.

C3

The simulated wave had to break through Cenotaph Island, opening a narrow channel through the trees (Miller, 1960).

C4

The trimline maximum distance of 1100 m from the high-tide shoreline at Fish Lake had to be reached.

More than 250 simulations were performed in order to find the optimal values for the parameters that best verified the four conditions mentioned above. Finally, the optimal parameters selected were

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathit{\alpha }=\mathrm{13}{}^{\circ },r=\mathrm{0.44},m_{\mathrm{f}}=\mathrm{0.08}.\end{array}$$

Setting the three parameters to the values given above, simulation satisfied the previous four criteria with very good accuracy, more precisely,

SC1

the run-up in Gilbert Inlet reached 523.85 m (Fig. 10a),

SC2

the run-up peak in the vicinity of Mudslide Creek reached a height close to 200 m (Fig. 10b),

SC3

the wave produced a narrow channel crossing through Cenotaph Island (Fig. 10c), and, finally,

SC4

the run-up reached more than 1100 m in distance from the high-tide shoreline in the Fish Lake area (Fig. 10c).

Figure 10 graphically shows the four selected criteria enumerated above. In addition, in Fig. 7 the maximum run-up values for the optimal solution in the four regions selected for the sensitivity analysis are marked in the four panels and can be compared with the observed reference values. The next section describes in some detail the numerical experiment performed with the optimal set of parameters.

Figure 12 Wave evolution from t=40 s until t=2 min. t=30 s. In panels (b, d, f, h) the track of the inundated areas is kept. (a, c, e, g) South view; (b, d, f, h) plan view.

9 Model results

Sensitivity analysis provided the optimal set for the three key parameters considered for this study. In this section, model results corresponding to the simulation performed with these optimal parameters are presented: first the main characteristics of the giant wave generated in Gilbert Inlet, second optimal description of wave evolution through the main stem of the Lituya Bay, and third inundation details resulting from a comparison of numerical simulation run-up with the real trimline observed in several areas of interest. Additional material in the form of a numerical simulation movie is also provided.

Figure 13 Wave evolution from t=3 min and 30 s until t=7 min. In panels (b, d, f, h) the track of the inundated areas is kept. (a, c, e, g) South view; (b, d, f, h) plan view.

9.1 Giant wave generation and evolution in Gilbert Inlet, t: 0–39 s

Following the landslide trigger, the generated wave reaches its maximum amplitude, 272.4 m, at t=8 s (Fig. 11b). The wave spreads outwards in a southward direction with decreasing amplitude (Fig. 11c, d). A total of 22 s after triggering, the wave hits the bottom of the southwest spur of Gilbert Inlet. At t=39 s the maximum 523.9 m run-up is reached (Fig. 11f).

9.2 Wave evolution

9.2.1 t: 30 s–2 min

While the maximum run-up on the east side of Gilbert head is reached, the southern propagating part of the initial wave, with a height of more than 100 m, moves in a southwest direction, hitting the south shoreline of the bay after 35 s (Fig. 12a, b). The impact causes maximum run-up close to 180 m to occur in the vicinity of Mudslide Creek at t=70 s (Fig. 12c, d). In the meanwhile, part of the water reaching the maximum run-up area over Gilbert head retreats and part flows over Gilbert head, inundating the observed affected area to the south.

Figure 14 Simulated inundated areas and maximum height all around Lituya Bay. Trimline in pink.

9.2.2 t: 2–3 min

While the initial wave propagates through the main axis of Lituya Bay to Cenotaph Island, a larger second wave appears as a reflection of the first one from the south shoreline (Fig. 12c, d). Both waves swept each of the shorelines in their path. Along the north shoreline, wave run-up reaches between 50 and 80 m in height (Fig. 12e), while along the south shoreline the run-up reaches heights between 60 and 150 m.

The first wave reaches Cenotaph Island after 2 min and 5 s with a mean amplitude of close to 20 m (Fig. 12e, f), flooding over more than 650 m from the most eastern prominence of the island and about 700 m from the little cape, slightly southwards. About 25 s later, a second wave of approximately 32 m in height hits the east coast (Fig. 12g, h).

9.2.3 t: 3–5 min

After hitting Cenotaph Island, the wave splits into two parts; one advances in the shallow channel north of the island and the second travels through the deeper channel south of the island (Fig. 13a, b). Waves higher than 25 m hit the north shoreline area in front of Cenotaph Island causing large extent run-ups (above 1 km inland from the coastline) in the east area near Fish Lake (Fig. 13c, d). Along the south shoreline in front of Cenotaph Island, larger waves with 40–50 m of amplitude hit the coast and penetrate about 1 km inland on the flat areas located east of The Paps (Fig. 13e, f).

9.2.4 t: 5–8 min

Large inundated areas were formed both around Fish Lake, on the north shoreline, and in the flat areas surrounding The Paps, while the main wave reaches the narrow area near La Chaussee Spit. During this time, the wave amplitude is larger than 15 m, indeed, over 20 m on the north shoreline. At 5 min and 50 s the wave reaches La Chaussee Spit, passing over before reaching the sea and then partially reflecting a wave back into the bay (Fig. 13g, h).

9.3 Inundation assessment

In this section, the computed inundated area limit is compared with the real trimline drawn by Miller. Figure 14 depicts the total extent of the inundated area and the maximum water height reached. The trimline determined by Miller (Fig. 1) is superimposed in pink. In order to closely assess the accuracy of the predicted numerical inundated area with respect to the observed trimline, we have performed this comparison focusing successively on the different areas of interest.

9.3.1 Gilbert Inlet

As it has been shown, on Gilbert Head, the maximum run-up (523.85 m) is reached on the east slope. Furthermore, the run-up is extended oblique to slope on the western face of Gilbert Head. Run-up extent and trimline coincide quite remarkably in this sector (Fig. 10a). There is large inundated area over the Lituya Glacier, from its shoreline up to more than 2 km over the glacier. Finally, a good correlation between trimline and inundation on the east slope of Gilbert Inlet where the slide was initially located is observed.

9.3.2 North shoreline

Figure 14 shows the good agreement between the model-simulated inundated area and the real trimline around the east part of the north shore due to the higher slopes. Good agreement between trimline and simulated run-up on the north shoreline of La Chaussee Spit is also found.

9.3.3 Fish Lake area

The agreement between model and observation around the flat areas surrounding Fish Lake is good (Fig. 10c). Here, the inundation extent includes vicinity areas of Fish Lake under 40 m in height. In order to achieve a better agreement between simulated inundation limit and observed trimline, it would be necessary to consider a map of drag friction to account for different friction coefficients depending on the type of vegetation or soil (Kaiser et al., 2011). In flat areas, good agreement between model results and observed data requires that the presence of vegetation or any other type of obstacle not collected in the topography be taken into account or parameterized in some way.

9.3.4 South shoreline

The computed run-up underestimates the trimline over steeper slopes on the eastern third of the south shores (Fig. 14). Moreover, the numerical model provides mean run-up heights of around 120 m while trimline heights move from 140 to 200 m around this sector. The reason for this mismatch is probably due to numerical resolution. In order to capture the steep slopes in this area, a higher numerical resolution than provided by available data would be required. There is good agreement along The Paps shores as noted on the south shoreline in front of La Chaussee Spit (Fig. 10d).

9.3.5 East area of The Paps

Around the central third of the south shoreline, the agreement between model and observed trimline is less coherent (Fig. 10d). On the flat east area near The Paps, large inundation occurs, flooding the sector area with a height of 30–40 m. As already mentioned, it would be necessary to consider a map of drag frictions in order to achieve more precise results over this area.

9.3.6 Cenotaph Island

One of the items to be checked in the sensitivity analysis presented in Sect. 8 was related to flooding on Cenotaph Island. Therefore, good agreement is expected in this particular location. In fact, the trimline and the computed inundation limit closely track on the island as shown in Fig. 10c.

9.3.7 La Chaussee Spit

As has been previously described, La Chaussee Spit was completely covered by water from minute 6 for more than 90 s (Figs. 13g, h). The computed inundation around La Chaussee Spit is in good agreement with the observed trimline (Fig. 10d).

Figure 15 Same as in Fig. 10 but for the blind simulation with standard not optimally adjusted parameters. Comparison with the observed trimline. (a) Gilbert and Crillon inlets. (b) Inner part of Lituya Bay from Cenotaph Island. (c) Fish Lake and Cenotaph Island area. (d) Entrance of the bay up to Cenotaph Island.

10 Discussion

10.1 Potential sources of error

The Landslide-HySEA model was tested against analytic solutions and laboratory measurements of Fritz et al. (2009) (Sánchez-Linares, 2011). In addition, Landslide-HySEA recently participated in the “NTHMP/MMS Landslide Model Benchmarking Workshop” hosted by Texas A&M University at Galveston on 9–11 January 2017 during which the model satisfactorily reproduced expected results. Landslide-HySEA is a finite-volume coupled landslide–fluid model that acts as a shallow water model when the slide layer is immobile or when there is not a sediment layer in the column of water. A suitable and simple treatment of the wet–dry fronts avoids spurious wave reflection on the coast and produces a realistic simulation of flooded areas. Transitions between sub- and supercritical flows, that continuously appear in simulations including those presented here and which further complicate matters, are also suitably treated.

Nonetheless, as it was previously mentioned, potential sources of model errors are the quality of model initialization parameters, the initial landslide conditions, or the digital elevation model (DEM) due to limitations associated with bathymetric data. Moreover, in real landslides, the material is neither homogeneous nor granular, as assumed here in the present study. Nevertheless, this type of model can be used in practice to provide general information about the generated tsunami and the flooded areas as demonstrated in the presented results.

10.1.1 Limitations of the DEM and digital bathymetry

A high-quality DEM is necessary to properly model tsunami wave dynamics and inundation onshore, especially in areas with complicated bathymetries. In this study there was an additional need for good information on topo-bathymetric data just before the event in order to produce realistic pre-event geometry.

Though we have combined the DEM based on the best available data in the region (described in Sect. 4), neither pre-tsunami bathymetry data of the bay nor the definition of Lituya Glacier front just before the 1958 tsunami were available with fine enough resolution or quality. Of even greater significance, estimations of the volume and position of the slide that caused the tsunami were all that were available. Thus, as was described in Sects. 3 and 5, a proposed reconstruction of the original Lituya Glacier shoreline provided by Miller (1960) and data from the 1926 and 1959 US Coast and Geodetic Surveys were used in this study.

10.2 Model results

Due to the choice of optimal parameters in the sense described in Sect. 8, the simulation presented achieves the main objectives proposed. Section 9 presents the first stage of the tsunami dynamics in Gilbert Inlet: giant wave generation and the inundation induced over the east slope of Gilbert Inlet. Later, the wave propagation in the southwest direction along the Lituya Bay is described until the wave crosses La Chaussee Spit. Model results are in good agreement with those described in Miller (1960).

In a second study stage, an inundation assessment is performed. A detailed description of the run-up areas along the shores of the bay is presented. In general, computed inundation areas are in very good agreement with Miller observations. Nevertheless, the model provides larger inundation areas than the 10.35 km² between the trimlines and the high-tide shorelines estimated by Miller. It is noted, however, that Miller made an estimate of the total area inundated by the wave of at least 13 km², an estimate that is closer to the model results.

11 Conclusions and future work

It has been demonstrated that the landslide-triggering mechanism proposed by Fritz et al. (2009) is crucial in order to reproduce not only the wave dynamics inside Gilbert Inlet, but also all tsunami dynamics produced along the bay, including inundation effects, wave heights, and several details observed in Miller (1960). The simulated wave heights and run-up (as assessed by the trimline location) are in good agreement with the majority of observations and conclusions described by Miller (1960).

It has been shown that the numerical model used can simulate subaerial scenarios similar to the Lituya Bay case provided that some information is available to calibrate the model. The main question that remains to be answered is then obvious: what happens when information to calibrate the model is not available? In that case, which approach is followed? In other words, how would an actual risk assessment study would be performed without post-event information? In that case two approaches can be followed. One first option is a deterministic approach, in which, depending on the characteristics of the slide, some coefficients are selected as default. In the case considered in this work, consisting in a slide moving, essentially as a solid block, the “blind” proposed parameters would be (in parentheses the optimal values for comparison)

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathit{\alpha }=\mathrm{16}{}^{\circ }\left(\mathrm{13}{}^{\circ }\right),r=\mathrm{0.3}\left(\mathrm{0.44}\right),m_{\mathrm{f}}=\mathrm{0.08}\left(\mathrm{0.08}\right).\end{array}$$

The numerical results for this case, which must be compared with the optimal case selected in Fig. 10, are presented in Fig. 15. As can be observed, in this case, blind numerical results produce higher run-up and larger inundated areas than the optimal calibrated simulation; therefore, from the point of view of risk assessment, considering the standard parameters should not underestimate the risk. Furthermore, in some areas the agreement is better than for the optimal case, although the four required criteria at the same time are better achieved in the optimal case.

Concerning future work, as uncertainty in the data (initial condition, model parameters, etc.) is of paramount importance in real applications, a promising line of research is uncertainty quantification. Therefore, some information of the main probabilistic moments should be provided. Uncertainty quantification is currently a very active area of research, with one of the most efficient techniques utilized being multilevel Monte Carlo methods. To run such a method, a family of embedded meshes is first considered. Then, a large enough number of samples of the stochastic terms are chosen and, for each sample, a deterministic simulation is run. Finally, the probabilistic moments are then computed by a weighted average of the deterministic computations (Sánchez-Linares et al., 2016).

Another improvement of the model envisioned will be carried out by considering shallow Bingham dense avalanche models, like those introduced in Fernández-Nieto et al. (2010), that will be coupled with the hydrodynamic model. Conversely, current efforts are focused on the implementation of a model including dispersive effects to allow for comparison of model performance with and without dispersion (Macías et al., 2017a). This effort too will serve to assess the role played by dispersive terms in these landslide-generated types of events. In any case, the present work shows that a Savage–Hutter model coupled with shallow water equations is sufficient to suitably reproduce the main features of an extreme event such as the one occurring in Lituya Bay in 1958.