Research article 15 Feb 2019
Research article  15 Feb 2019
The Lituya Bay landslidegenerated megatsunami – numerical simulation and sensitivity analysis
 ^{1}Dpto. de Matemática Aplicada, ETSII, Universidad de Málaga, 29080, Málaga, Spain
 ^{2}Dpto. de Análisis Matemático, Estadístice e Invetigación Operativa y Matemática Aplicada, Facultad de Ciencias, Universidad de Málaga, 29080, Málaga, Spain
 ^{3}Unit of Numerical Methods, SCAI, Universidad de Málaga, 29080, Málaga, Spain
 ^{4}NOAA/Pacific Marine Environmental Laboratory (PMEL), Seattle, WA, USA
 ^{1}Dpto. de Matemática Aplicada, ETSII, Universidad de Málaga, 29080, Málaga, Spain
 ^{2}Dpto. de Análisis Matemático, Estadístice e Invetigación Operativa y Matemática Aplicada, Facultad de Ciencias, Universidad de Málaga, 29080, Málaga, Spain
 ^{3}Unit of Numerical Methods, SCAI, Universidad de Málaga, 29080, Málaga, Spain
 ^{4}NOAA/Pacific Marine Environmental Laboratory (PMEL), Seattle, WA, USA
Correspondence: Jorge Macías (jmacias@uma.es)
Hide author detailsCorrespondence: Jorge Macías (jmacias@uma.es)
The 1958 Lituya Bay landslidegenerated megatsunami is simulated using the LandslideHySEA model, a recently developed finitevolume Savage–Hutter shallow water coupled numerical model. Two factors are crucial if the main objective of the numerical simulation is to reproduce the maximal runup with an accurate simulation of the inundated area and a precise recreation of the known trimline of the 1958 megatsunami 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 topobathymetric data, and on information extracted from the work of Miller (1960), an approximation of Gilbert Inlet topobathymetry was set up and used for the numerical simulation of the megaevent. 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 megatsunami. 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 3D basin geometry, where no smoothing or other stabilizing factors in the bathymetric data are applied.
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 socalled megatsunamis, which are characterized by localized extreme runup 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 nondiffusive linear approximation. With landslidegenerated 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 subaeriallandslidegenerated 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 runup 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 runup 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 runup 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 runup produced by an aerial landslide in an enclosed bay, a simulation that has not been successfully undertaken previously with more comprehensive numerical models.
Full 3D numerical modeling of landslidegenerated tsunamis (Horrillo et al., 2013) is uncommonly used for realworld scenarios due to the highly demanding computational resources required. Whereas common thought persists that NLSW equations are sufficient for simulation of oceanwide tsunami propagation averaged models, the importance of frequency dispersion for modeling landslidegenerated 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 landslidegenerated 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 landslidegenerated tsunamis. For the case of enclosed basins, bays, and fjords, it was agreed that NLSW models remain a suitable tool.
Among all examples of subaeriallandslidegenerated 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 2D 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ánchezLinares, 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 runup 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 megatsunami in a realistic threedimensional geometry with no smoothing in the geometry or initial conditions.
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^{3} in Gilbert Inlet (Miller, 1960) that in turn produced the largest tsunami runup 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 PararasCarayannis (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 twodimensional scaled physical model of the Gilbert Inlet. A pneumatic landslide generator was used to generate a highspeed granular slide with density and volume based on Miller (1960) impacting the water surface at a mean velocity of 110 m s^{−1}. The experimental runup matches the trimline of forest destruction on the spur ridge in Gilbert Inlet.
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^{−1} slide impact velocity that Fritz et al. (2009) measure in their experiment. In the same way that the laboratory experiment reproduced the observed runup, it is the way in which the numerical experiments performed here have been initialized.
Lituya Bay (Fig. 1), located within Glacier Bay National Park on the northeast shore of the Gulf of Alaska, is a Tshaped 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, fjordlike, 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.
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).
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 Ushaped 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 preevent 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 pretsunami 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^{6} m^{3}, that is, very close to the slide volume estimated by Miller (1960).
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^{6} m^{3}.
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^{6} m^{3}, 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^{−1}, which is the approximated impact velocity between the slide and the water surface estimated by Fritz et al., assuming freefall 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^{−1} 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^{−1}. 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 runup 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^{−1}.
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^{−1} (see Fig. 3). Thus, with this initial condition, the model reproduces both the runup 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.
Coulombtype models for granulardriven flows have been intensively investigated in the last decade, following the pioneering work of Savage and Hutter (Savage and Hutter, 1989), who derived a shallowwatertype 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 finitevolume numerical model, the LandslideHySEA 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 twolayer Savage–Hutter model introduced in FernándezNieto 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ánchezLinares (2011), the LandslideHySEA model is used to reproduce the Fritz et al. (2009) laboratory experiment, producing good results.
This section describes the system of partial differential equations modeling landslidegenerated tsunamis based on layered average models. The 2D LandslideHySEA model is a twodimensional version of the model proposed in FernándezNieto et al. (2008) for 1D problems, for which local coordinates are not considered.
6.1 Simplified twolayer Savage–Huttertype model
Let us consider a layered medium composed by a layer of inviscid fluid with constant homogeneous density ρ_{1} (water) and a layer of granular material with density ρ_{s} and porosity ψ_{0}. 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}})\phantom{\rule{0.25em}{0ex}}{\mathit{\rho}}_{\mathrm{s}}+{\mathit{\psi}}_{\mathrm{0}}{\mathit{\rho}}_{\mathrm{1}}$. The system of partial differential equations (PDEs) describing the coupled twolayer system is written as
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)\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.125em}{0ex}}\in \phantom{\rule{0.125em}{0ex}}D\phantom{\rule{0.125em}{0ex}}\subset \phantom{\rule{0.125em}{0ex}}{\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 nonerodible 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),\phantom{\rule{0.125em}{0ex}}{q}_{i,y}(x,y,t)\phantom{\rule{0.25em}{0ex}}\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)\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\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 (ρ_{1}<ρ_{2}). Note that H(x,y) does not depend on t, that is, the nonerodible 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_{1}, h_{2}, 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},\phantom{\rule{0.125em}{0ex}}{\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
The term ${\mathit{S}}_{c}\left(\mathit{W}\right)=\left({S}_{{c}_{x}}\right(\mathit{W}),\phantom{\rule{0.125em}{0ex}}{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
where m_{f} is a positive constant, and ${\mathit{S}}_{a}\left(\mathit{W}\right)=\left({S}_{{a}_{x}}\right(\mathit{W}),\phantom{\rule{0.125em}{0ex}}{S}_{{a}_{y}}(\mathit{W}\left)\right)$ parameterizes the friction between the fluid and the nonerodible bottom and is given by a Manning law (Dyakonova and Khoperskov, 2018):
where n_{1}>0 is the Manning coefficient.
${\mathit{S}}_{\mathrm{b}}\left(\mathit{W}\right)=\left({S}_{{\mathrm{b}}_{x}}\right(\mathit{W}),\phantom{\rule{0.125em}{0ex}}{S}_{{\mathrm{b}}_{y}}(\mathit{W}\left)\right)$ parameterizes the friction between the granular and the nonerodible bottom and as in the previous case is given by a Manning law:
where n_{2}>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_{2}=0. Similarly, if ${h}_{\mathrm{1}}(x,y,t)=\mathrm{0}$ then m_{f}=0 and n_{1}=0.
Finally, $\mathit{\tau}=({\mathit{\tau}}_{x},\phantom{\rule{0.125em}{0ex}}{\mathit{\tau}}_{y})$ is defined as follows:
where ${\mathit{\sigma}}^{c}=g(\mathrm{1}r)\phantom{\rule{0.25em}{0ex}}{h}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{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ándezNieto et al., 2008). In the next section, the finitevolume scheme used to discretize system Eq. (1) is described. As friction terms are semiimplicitly 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.
To discretize system Eq. (1), the domain D is divided into L cells or finite volumes V_{i}⊂ℝ^{2}, $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}∈ℝ^{2} 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},\phantom{\rule{0.125em}{0ex}}{\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^{n}:
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 onedimensional 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ándezNieto 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:
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ánchezLinares 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)$, $\phantom{\rule{0.125em}{0ex}}{S}_{{\mathrm{f}}_{\mathrm{2}}}\left(\mathit{W}\right)$, $\phantom{\rule{0.125em}{0ex}}{S}_{{\mathrm{f}}_{\mathrm{3}}}\left(\mathit{W}\right)$, and ${S}_{{\mathrm{f}}_{\mathrm{4}}}\left(\mathit{W}\right)$ are discretized semiimplicitly; next the Coulomb friction term τ will be discretized following FernándezNieto et al. (2008). The resulting numerical scheme is a threestep 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 threestep method will be denoted as follows:
The resulting scheme is exactly well balanced (SánchezLinares et al., 2015) for the stationary water at rest solution (${\mathit{q}}_{\mathrm{1}}={\mathit{q}}_{\mathrm{2}}=\mathbf{0}$ and μ_{1} and μ_{2} constant). Moreover, the scheme accurately solves the stationary solutions corresponding to ${\mathit{q}}_{\mathrm{1}}={\mathit{q}}_{\mathrm{2}}=\mathbf{0}$, μ_{1} constant and ${\partial}_{x}{\mathit{\mu}}_{\mathrm{2}}<\mathrm{tan}\left(\mathit{\alpha}\right)$ and ${\partial}_{y}{\mathit{\mu}}_{\mathrm{2}}<\mathrm{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 GPUbased 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.
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:
The other two parameters required in model parameterizations, the Manning coefficients n_{1} and n_{2}, were set constant with standard values of n_{1}=0.02 and n_{2}=0.05 (Arcement and Schneider, 1989; Phillips and Tadayon, 2006).
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 runup 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 3by4 graphs presented here as three figures (Figs. 7, 8, and 9).
Each panel in Fig. 7 summarizes the variability in the maximum computed runup 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 runup for different values of the friction coefficient, m_{f}, is presented. Finally, the vertical bars show the runup 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 runup 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 runup 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 runup 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 runups for different values of the Coulomb friction angle, α. Finally, the vertical bars show the runup 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.
In Fig. 9, each panel shows the values of the maximum runup 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 runups for different values of the friction coefficient, m_{f}. Finally, the vertical bars show the runup 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.
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.
Then, the four criteria considered in order to select the optimal parameters were as follows.
 C1

The runup 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 hightide 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
Setting the three parameters to the values given above, simulation satisfied the previous four criteria with very good accuracy, more precisely,
 SC1

the runup in Gilbert Inlet reached 523.85 m (Fig. 10a),
 SC2

the runup 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 runup reached more than 1100 m in distance from the hightide 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 runup 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.
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 runup with the real trimline observed in several areas of interest. Additional material in the form of a numerical simulation movie is also provided.
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 runup is reached (Fig. 11f).
9.2 Wave evolution
9.2.1 t: 30 s–2 min
While the maximum runup 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 runup 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 runup area over Gilbert head retreats and part flows over Gilbert head, inundating the observed affected area to the south.
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 runup reaches between 50 and 80 m in height (Fig. 12e), while along the south shoreline the runup 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 runups (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 runup (523.85 m) is reached on the east slope. Furthermore, the runup is extended oblique to slope on the western face of Gilbert Head. Runup 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 modelsimulated 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 runup 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 runup underestimates the trimline over steeper slopes on the eastern third of the south shores (Fig. 14). Moreover, the numerical model provides mean runup 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.
10.1 Potential sources of error
The LandslideHySEA model was tested against analytic solutions and laboratory measurements of Fritz et al. (2009) (SánchezLinares, 2011). In addition, LandslideHySEA 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. LandslideHySEA is a finitevolume 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 highquality 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 topobathymetric data just before the event in order to produce realistic preevent geometry.
Though we have combined the DEM based on the best available data in the region (described in Sect. 4), neither pretsunami 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 runup 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^{2} between the trimlines and the hightide 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^{2}, an estimate that is closer to the model results.
It has been demonstrated that the landslidetriggering 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 runup (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 postevent 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)
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 runup 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ánchezLinares 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ándezNieto 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 landslidegenerated 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.
Our definition of the bathymetry is based on public domain topobathymetric data cited in the paper and on information extracted from the work of Miller (1960). Data concerning the numerical results can be obtained by a personal request to the corresponding author. Any material requested from the authors will be made available in a public repository of our university. A visualization of the numerical simulation is available as additional material at https://doi.org/10.5446/39493 (Macías et al., 2019c).
JMGV designed, with available topobathymetric data and literature, the initial condition and wrote the first version of the paper. JM wrote the subsequent versions of the paper and carried out the entire revision process, performed the sensitivity analysis, and generated the corresponding figures. MJC contributed to the definition of the mathematical model and the design of the numerical scheme. CSL generated the reconstructed topobathymetry used in the simulations and developed the discretization of the source terms. SOA performed the numerical simulations and produced some of the figures. MA implemented the CUDA GPU version of the numerical code used. He also performed some numerical tests and provided some figures. DA collected the historical information officially available for this event in the USGS and NOAA archives. He contributed to the writing of the initial version of the article.
The authors declare that they have no conflict of interest.
All numerical experiments required for the development of this research have been performed at the
Unit of Numerical Methods of the University of Málaga (Ada Byron Building).
This work was partially funded by the NOAA Center for Tsunami Research (NCTR), Pacific Marine Environmental
Laboratory (USA) contract no. WE133R12SE0035, by the Junta de Andalucía research project TESELA (P11RNM7069),
by the Spanish Government Research project SIMURISK
(MTM201570490C0201R), and by the Universidad de Málaga, Campus de Excelencia
Andalucía TECH. All data required to perform the numerical simulations
presented in this study are described in the text or are publicly available.
Diego Arcas was supported by the National Oceanic and Atmospheric Administration
with this report being PMEL contribution no. 4665.
Edited by: Maria Ana Baptista
Reviewed by: two
anonymous referees
Abadie, S., Gandon, C., Grilli, S., Fabre, R., Riss, J., Tric, E., and Morichon, D.: 3D numerical simulations of waves generated by subaerial mass failures. Application to La Palma Case, in: Proc. 31st Intl. Coastal Engng. Conf., edited by: Mc Kee Smith, J., ICCE08, Hamburg, Germany, September, 2008, 1384–1395, World Scientific Publishing Co. Pte. Ltd., https://doi.org/10.1142/9789814277426_0115, 2009. a
Arcement, G. J. and Schneider, V. R.: Guide for Selecting Manning's Roughness Coefficients for Natural Channels and Flood Plains, USGS WaterSupply papers, 2339, 44 pp., 1989. a
Basu, D., Das, K., Green, S., Janetzke, R., and Stamatakos, J.: Numerical simulation of surface waves generated by subaerial landslide at Lituya Bay Alaska, J. Offshore Mech. Eng., 132, 11, https://doi.org/10.1115/1.4001442, 2010. a
Bloom, C. K., Macinnes, B., Higman, B., Stark, C. P., Lynett, P., Ekström, G., Hibert, C., Willis, M. J., and Shugar, D. H.: Field observations from a massive landslide tsunami in Taan Fjord, Wrangell St. Elias National Park, AK, Geological Society of America Abstracts with Programs, vol. 48, https://doi.org/10.1130/abs/2016AM285078, GSA Annual Meeting in Denver, Colorado, USA, 2016. a
Bouchut, F. and Westdickenberg, M.: Gravity driven shallow water model for arbitrary topography, Comm. Math. Sci., 2, 359–389, 2004. a
Bouchut, F., Mangeney, A., Perthame, B., and Vilotte, J. P.: A new model of Saint Venant and SavageHutter type for gravity driven shallow flows, C. R. Acad. Sci. Paris, Ser I, 336, 531–536, 2003. a
Bouchut, F., FernándezNieto, E. D., Mangeney, A., and Lagrée, P.Y.: On new erosion models of Savage–Hutter type for avalanches, Acta Mech., 199, 181–208, https://doi.org/10.1007/s0070700705349, 2008. a
Bouchut, F., FernándezNieto, E. D., Mangeney, A., and NarbonaReina, G.: A twophase twolayer model for fluidized granular flows with dilatancy effects, J. Fluid Mech., 801, 166–221, 2016. a
Castro, M. J., Ferreiro, A., GarcíaRodríguez, J. A., GonzálezVida, J. M., Macías, J., Parés, C., and VázquezCendón, M. E.: The numerical treatment of wet/dry fronts in shallow flows: application to onelayer and twolayer systems, Math. Comput. Model., 42, 419–439, 2005. a, b
Castro, M. J., Gallardo, J. M., and Parés, C.: High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products, Applications to shallowwater systems, Math. Comp., 75, 1103–1134, 2006. a
Castro, M. J., ChacónRebollo, T., FernándezNieto, E., GonzálezVida, J. M., and Parés, C.: Wellbalanced finite volume schemes for 2d nonhomogeneous hyperbolic systems, Application to the dam break of Aznalcóllar, Comput. Meth. Appl. Mech. Eng., 197, 3932–3950, 2008. a
Castro, M. J., FernándezNieto, E., Morales, T., NarbonaReina, G., and Parés, C.: A HLLC scheme for nonconservative hyperbolic problems, Application to turbidity currents with sediment transport, ESAIM: Math. Model. Numer. Anal., 47, 1–32, 2012. a
de la Asunción, M., Castro, M. J., FernándezNieto, E., Mantas, J., Ortega, S., and GonzálezVida, J. M.: Efficient GPU implementation of a two waves TVDWAF method for the twodimensional one layer shallow water system on structured meshes, Comp. Fluids, 80, 441–452, 2013. a
de la Asunción, M., Castro, M. J., Mantas, J., and Ortega, S.: Numerical simulation of tsunamis generated by landslides on multiple GPUs, Adv. Eng. Softw., 99, 59–72, 2016. a
Dyakonova, T. A. and Khoperskov, A. V.: Bottom friction models for shallow water equations: Manning's roughness coefficient and smallscale bottom heterogeneity, J. Phys. Conf. Ser., 973, 1–10, 012032, https://doi.org/10.1088/17426596/973/1/012032, 2018. a
FernándezNieto, E., Bouchut, F., Bresch, D., Castro, M. J., and Mangeney, A.: A new Savage Hutter type model for submarine avalanches and generated tsunami, J. Comput. Phys., 227, 7720–7754, https://doi.org/10.1016/j.jcp.2008.04.039, 2008. a, b, c, d
FernándezNieto, E., Noble, P., and Vila, J.P.: Shallow water equations for nonnewtonian fluids, J. NonNewt. Fluid Mech., 165, 712–732, https://doi.org/10.1016/j.jnnfm.2010.03.008, 2010. a
FernándezNieto, E., Castro, M. J., and Parés, C.: On an Intermediate Field Capturing Riemann solver based on a Parabolic Viscosity Matrix for the twolayer shallow water system, J. Sci. Comput., 48, 117–140, 2011. a
Fritz, D., Hager, W., and Minor, H.E.: Lituya Bay case: rockslide impact and wave runup, Science of Tsunami Hazards, 19, 3–22, 2001. a
Fritz, H., Mohammed, F., and Yoo, J.: Lituya Bay landslide impact generated megatsunami 50th anniversary, Pure Appl. Geophys., 166, 153–175, 2009. a, b, c, d, e, f, g, h, i
Gallardo, J., Parés, C., and Castro, M. J.: On a wellbalanced highorder finite volume scheme for shallow water equations with topography and dry areas, J. Comput. Phys., 227, 574–601, 2007. a
Gray, J. M. N. T., Wieland, M., and Hutter, K.: Gravitydriven free surface flow of granular avalanches over complex basal topography, P. R. Soc. Lond. A, 455, 1841–1874, https://doi.org/10.1098/rspa.1999.0383, 1999. a
Harbitz, C., Pedersen, G., and Gjevik, B.: Numerical simulation of large water waves due to landslides, J. Hydraul. Eng., 119, 1325–1342, 1993. a
Higman, B., Shugar, D. H., Stark, C. P., Ekström, G., Koppes, M. N., Lynett, P., Dufresne, A., Haeussler, P. J., Geertsema, M., Gulick, S., Mattox, A., Venditti, J. G., Walton, M. A. L., McCall, N., Mckittrick, E., MacInnes, B., Bilderback, E. L., Tang, H., Willis, M. J., Richmond, B., Reece, R. S., Larsen, C., Olson, B., Capra, J., Ayca, A., Bloom, C., Williams, H., Bonno, D., Weiss, R., Keen, A., Skanavis, V., and Loso, M.: The 2015 landslide and tsunami in Taan Fiord, Alaska, Sci. Rep. UK, 8, 12993, https://doi.org/10.1038/s4159801830475w, 2018. a
Horrillo, J., Wood, A., Kim, G.B., and Parambath, A.: A simplified 3D NavierStokes numerical model for landslidetsunami: Application to the Gulf of Mexico, J. Geophys. Res.Oceans, 118, 6934–6950, https://doi.org/10.1002/2012JC008689, 2013. a
Jørstad, F.: Waves generated by slides in norwegian fjords and lakes, Tech. rep., Norwegian Geotechnical Institute, 1968. a
Kaiser, G., Scheele, L., Kortenhaus, A., Løvholt, F., Römer, H., and Leschka, S.: The influence of land cover roughness on the results of high resolution tsunami inundation modeling, Nat. Hazards Earth Syst. Sci., 11, 2521–2540, https://doi.org/10.5194/nhess1125212011, 2011. a
Lindstrøm, E., Pedersen, G., Jensen, A., and Glimsdal, S.: Experiments on slide generated waves in a 1:500 scale fjord model, Coast. Eng., 92, 12–23, https://doi.org/10.1016/j.coastaleng.2014.06.010, 2014. a
Løvholt, F., Pedersen, G., and Gisler, G.: Oceanic propagation of a potential tsunami from the La Palma Island, J. Geophys. Res., 113, C09026, https://doi.org/10.1029/2007JC004603, 2008. a
Løvholt, F., Glimsdal, S., Lynett, P., and Pedersen, G.: Simulating tsunami propagation in fjords with longwave models, Nat. Hazards Earth Syst. Sci., 15, 657–669, https://doi.org/10.5194/nhess156572015, 2015. a, b, c
Macías, J., FernándezSalas, L. M., GonzálezVida, J. M., Vázquez, J. T., Castro, M. J., Bárcenas, P., Díaz del Río, V., Morales de Luna, T., de la Asunción, M., and Parés, C.: Deslizamientos Submarinos y Tsunamis en el Mar de Alborán. Un Ejemplo de Modelización, Temas de Oceanografía, vol. 6, Instituto Español de Oceanografía, Madrid, 2012. a
Macías, J., Vázquez, J. T., FernándezSalas, L. M., GonzálezVida, J. M., Bárcenas, P., Castro, M. J., Díaz del Río, V., and Alonso, B.: The AlBorani submarine landslide and associated tsunami, A modelling approach, Mar. Geo., 361, 79–95, https://doi.org/10.1016/j.margeo.2014.12.006, 2015. a
Macías, J., Castro, M. J., Escalante, C., Ortega, S., and GonzálezVida, J. M.: HySEA model, NTHMP Landslide Benchmarking Results, Tech. rep., MMS/NTHMP, https://doi.org/10.13140/RG.2.2.27081.60002, 2017a. a, b
Macías, J., Castro, M. J., Ortega, S., Escalante, C., and GonzálezVida, J. M.: Performance benchmarking of the TsunamiHySEA model for NTHMP's inundation mapping activities, Pure Appl. Geophys., 174, 3147–3183, https://doi.org/10.1007/s0002401715831, 2017b. a
Macías, J., Castro, M. J., and Escalante, C.: Performance assessment of TsunamiHySEA model for NTHMP tsunami currents benchmarking, Part I Laboratory data, Coast. Eng., in review, 2019a. a
Macías, S., Castro, M. J., GonzálezVida, J. M., and Ortega, S.: Performance assessment of TsunamiHySEA model for NTHMP tsunami currents benchmarking, Part II Field data, Coast. Eng., in review, 2019b. a
Macías, J., de la Asunción, M., Ortega, S., GonzálezVida, J. M., and Castro, M. J.: Simulation of the 1958 Lituya Bay megatsunami, https://doi.org/10.5446/39493, 2019c.
Mader, C. L.: Modeling the 1958 Lituya Bay megatsunami, Science of Tsunami Hazards, 17, 57–68, 1999. a
Mader, C. L. and Gittings, M. L.: Modeling the 1958 Lituya Bay mega tsunami, II, Science of Tsunami Hazards, 20, 241–250, 2002. a
Miller, D.: Giant Waves in Lituya Bay, Alaska: A Timely Account of the Nature and Possible Causes of Certain Giant Waves, with Eyewitness Reports of Their Destructive Capacity, Professional paper, US Government Printing Office, 1960. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x
Montagna, F., Bellotti, G., and Di Risio, M.: 3D numerical modeling of landslidegenerated tsunamis around a conical island, Nat. Hazards, 58, 591–608, https://doi.org/10.1007/s1106901096890, 2011. a
PararasCarayannis, G.: Analysis of mechanism of tsunami generation in Lituya Bay, Science of Tsunami Hazards, 17, 193–206, 1999. a
Pelanti, M., Bouchut, F., and Mangeney, A.: A RoeType Scheme for TwoPhase Shallow Granular Flows with Bottom Topography, Math. Model. Numer. Anal. (ESAIM:M2AN), 48, 851–885, https://doi.org/10.1051/m2an:2008029, 2008. a, b
Phillips, J. V. and Tadayon, S.: Selection of Manning's roughness coefficient for natural and constructed vegetated and nonvegetated channels, and vegetation maintenance plan guidelines for vegetated channels in central Arizona, US Geological Survey Scientific Investigations Report 2006/5108, 41 pp., 2006. a
Pitman, E. B. and Le, I.: A twofluid model for avalanche and debris flows, Philos. T. Roy. Soc. A, 363, 1573–1601, 2005. a
Quecedo, M., Pastor, M., and Herreros, M.: Numerical modelling of impulse wave generated by fast landslides, Int. J. Numer. Meth. Eng., 59, 1633–1656, https://doi.org/10.1002/nme.934, 2004. a
SánchezLinares, C.: Simulación numérica de tsunamis generados por avalanchas submarinas: aplicación al caso de Lituya Bay, Master's thesis, University of Málaga, Spain, 2011. a, b, c
SánchezLinares, C., de la Asunción, M., Castro, M. J., Mishra, S., and Sukys, J.: Multilevel Monte Carlo finite volume method for shallow water equations with uncertain parameters applied to landslidesgenerated tsunamis, Appl. Math. Model., 39, 7211–7226, 2015. a, b
SánchezLinares, C., Castro, M. J., de la Asunción, M., GonzálezVida, J. M., Macías, J., and Mishra, S.: Uncertainty quantification in tsunami modeling using Monte Carlo finite volume method, J. Math. Ind., 6, 26 pp., https://doi.org/10.1186/s1336201600228, 2016. a
Savage, S. and Hutter, K.: The motion of a finite mass of granular material down a rough incline, J. Fluid Mech., 199, 177–215, https://doi.org/10.1017/S0022112089000340, 1989. a, b
Sharpe, C.: Landslides and Related Phenomena, Columbia Univ. Press, New York, 1938. a
Schwaiger, H. F. and Higman, B.: Lagrangian hydrocode simulations of the 1958 Lituya Bay tsunamigenic rockslide, Geochem. Geophys. Geosyst., 8, Q07006, https://doi.org/10.1029/2007GC001584, 2007. a
US Coast and Geodesic Survey: Survey id: H04608: NOS Hydrographic Survey, 19261231, available at: https://data.world/usnoaagov/f6786b28ea064c9aac3053cb5356650c (last access: 7 February 2019), 1926. a
US Coast and Geodesic Survey: Survey id: H08492: NOS Hydrographic Survey, Lituya Bay, Alaska, 19590827, available at: https://data.world/usnoaagov/9401821a28f5484688db43e702a5b12b (last access: 7 February 2019), 1959. a
Varnes, D.: Landslides and Engineering Practice, Highhw. Res. Board Spec. Rep., vol. 29, chap. Landslides types and processes, 22–47, National Research Council (US), 1958. a
Weiss, R., Fritz, H. M., and Wünnemann, K.: Hybrid modeling of the megatsunami runup in Lituya Bay after half a century, Geophys. Res. Lett., 36, L09602, https://doi.org/10.1029/2009GL037814, 2009. a