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 re-creation of the known trimline of the 1958 mega-tsunami of Lituya Bay. First, the accurate reconstruction of the initial slide. Then, the choice of a suitable coupled landslide-fluid model able to reproduce how the energy released by the 5 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, this 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 10 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 3D basin geometry, where no smoothing or other stabilizing factors in the bathymetric data are applied.


Introduction
Tsunamis are most often generated by bottom displacements due to earthquakes.However, landslides, either submarine or subaerial, can also trigger devastating tsunami waves.Besides 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, October 17, 2015 (Bloom et al., 2016;Higman et al., 2017), or the very recent event of Sulawesi, Sept 29, 2018, among many others).
For seismic tsunami simulations, in general, the most critical phases are generation and arrival to a coast, including inundation.Propagation over deep basins can be modeled using the non-linear shallow water (NLSW) equations or more typically using a non-diffusive linear approximation.With landslide generated tsunamis, however, matters get 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 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 are 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 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 preferrence 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, on the other hand, is mainly responsible for coastal impact.(Lindstrøm et al., 2014) using a scaled laboratory set up 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 where 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 play against the argument of ruling out non-dispersive coupled models, such as the one proposed in the present work.Very recently, in February 2017 during the "2017 U.S. 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 Subcommitee (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 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á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 meter 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.

Background
At 06:16 UTC on July 10, 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 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 minutes according to two eyewitnesses who were anchored at the entrance of the bay.According to Miller (1960), between 1 and 2.5 minutes 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 by water of forest and vegetation (known as trimline) extended to a maximum of 524 m above mean sea level on the spur southwest of Gilbert Inlet (Figure 1).Maximum inundation distance reached to 1,400 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:1,000 scale was constructed at the University of California, Berkeley (R.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 three 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).Geologic Survey, 1961) showing the settings and trimline of 1958 megatsunami (based on data from Miller (1960)).Units in feet.Key locations as Gilbert and Crillon Inlets, Cenotaph Island, Fish Lake or The Paps are shown.
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) and Fritz et al. (2009), the authors proposed and, in fact show, a landslide typology that, based on the generalized Froude similarity, make 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 runup matches the trimline of forest destruction on the spur ridge in Gilbert Inlet.Based on the experimental work of Fritz et al, 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 get the 110 m/s slide impact velocity that Fritz et al 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 has been initialized.

Area of Study
Lituya Bay (Figure 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 width ranging from 1.2 to 3.3 km except at the entrance, where the width is approximately 300 m.The north-eastward 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 m to 1,040 m in the foothills inmediately to the north and south, and more than 1,800 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 m and 1,290 m, respectively.

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 Figure 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 re-create 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 Figure 2).

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 (U.S. 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 meters 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 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 (U.S. 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 U.S. 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 Figures 2.c and 2.d it is shown the 1926 reconstructed bathymetry and the original 1959 bathymetry, respectively.The mass volume difference between these two bathymetries in Gilbert Inlet area is about 31 ⇥ 10 6 m 3 , that is, very close to the slide volume estimated by Miller (1960).

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 m to 915 m along the slope, a maximum thickness of about 92 m normal to the slope and a center of gravity at about 5 610 m 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 high, defining the upper bound for the mass slide.The volume of the reconstructed slide was of 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.

Landslide setup
In order to reproduce the main features of the slide impact, H. Fritz and collaborators 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, that it 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 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.
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 Figure 2.c), with a larger deposit concentration to the south part of the inlet was observed.
The observation of this fact made us to consider a slide initial velocity vector, v s , slightly shifted to the south, with modulus close to 80 m/s (see Figure 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.

Model description
Coulomb-type models for granular driven flows have been intensively investigated in the last decade, following the pioneering work of Savage-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), Bouchut andWestdickenberg (2004), Pelanti et al (2008), Bouchut et al. (2008), Bouchut et al. (2016)).In this framework, EDANYA group has implemented a finite volume numerical model, the Landslide-HySEA model (Castro et al., 2005(Castro et al., , 2006;;Gallardo et al., 2007;Castro et al., 2008Castro et al., , 2012;;Macías et al., 2012;de la Asunción et al., 2013;Macías et al., 2015), 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., 2018a(Macías et al., , 2017d) 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.
This section describes the system of partial differential equations modelling landslide generated tsunamis based on layered average models.The 2D Landslide-HySEA model, is a two-dimensional version of the model proposed in Fernández-Nieto et al. (2008) for 1D problems, where local coordinates are not considered.

Simplified Two-layer Savage-Hutter type 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 ⇢ 2 = (1 0 )⇢ s + 0 ⇢ 1 .The system of PDE describing the coupled two-layer system writes as: In these equations, subscript 1 refers to fluid upper layer, and subscript 2 to the lower layer composed of the fluidized material.i-th layer thickness at point (x, y) 2 D ⇢ R 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 q i (x, y, t) = (q i,x (x, y, t), q i,y (x, y, t)) is the flow of the i-th layer at point (x, y) at time t, that are related to the mean velocity of each layer (u i (x, y, t)) by q i (x, y, t) = h i (x, y, t) u i (x, y, t), i = 1, 2. The value r = ⇢ 1 /⇢ 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 non-erodible bottom topography does not change through the simulation although bottom may change due to second layer movement.Figure 5 graphically shows the relationship between h 1 , h 2 and H. Usually, h 1 + h 2 = H at rest or they represent the mean sea level.5 The term S c (W ) = S cx (W ), S cy (W ) 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 S a (W ) = S ax (W ), S ay (W ) parameterizes the friction between the fluid and the non-10 erodible bottom, and is given by a Manning law (Dyakonova and Khoperskov, 2018): where n 1 > 0 is the Manning coefficient.
S b (W ) = S bx (W ), S by (W ) parameterizes the friction between the granular and the non-erodible 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 2 (x, y, t) = 0.In this case, m f = 0 and n 2 = 0. Similarly, if h 1 (x, y, t) = 0 then m f = 0 and n 1 = 0.
System (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 (1) is described.As friction terms are semi-implicitly discretized, we first consider that S F (W ) = 0.Then, the way those terms are discretized is briefly described.
7 Numerical scheme To discretize system (1), the domain D is divided into L cells or finite volumes V i ⇢ R 2 , i = 1, . . ., L, which are assumed to be closed polygons.We assume here that the cells are rectangles with edges parallels to Cartesian axes.Given a finite volume V i , N i ⇢ R 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 ) is the unit normal vector to the edge ij and pointing towards V j .
We denote by W n i an approximation of the solution average over the cell V i at time t n : where |V i | is the area of cell V i and t n = t n 1 + t, where t is the time step.
Let us suppose that W n i 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 IFCP numerical scheme (see Fernández-Nieto et al. (2011)).Finally, W n+1 i is computed by averaging these approximate solutions.
The resulting numerical scheme writes as follows: To check the precise definition for the numerical fluxes, F ij (W n i , W n j , H i , H j ), see Sánchez-Linares et al. ( 2015) or de la Asunción et al. (2016).

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, that appears continuously in simulations 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.Coastline becomes a moving boundary, computed by the numerical scheme.Depending on the impact wave characteristics or the water back-flow movement, the computational cells are filled with water or they run dry, respectively.
Consequently, no specific stabilization model technique is either required.

Friction terms discretization
In this section the numerical scheme when the friction terms S F (W ) are discretized is presented.First, the terms S f1 (W ), S f2 (W ), S f3 (W ) y S f4 (W ) 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, where 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: The resulting scheme is exactly well-balanced Sánchez-Linares et al. ( 2015) for the stationary water at rest solution (q 1 = q 2 = 0 and µ 1 and µ 2 constant).Moreover, the scheme solves accurately the stationary solutions corresponding to q 1 = q 2 = 0, µ 1 constant and @ x µ 2 < tan(↵) and @ y µ 2 < tan(↵), that is a stationary water at rest solution for which the Coulomb friction term balances the pressure term in the granular material.

Numerical resolution and GPU implemetation
Landslide initial conditions have been described in Section 5.1.Initially water is at rest, and according to Miller (1960), an initial level of 1.52 m have been set.The computational grid considered for the numerical simulation is a 4 m ⇥ 7.5 m rectangular mesh composed of 3,650 ⇥ 1,271, i.e. 4,639,150, cells.The numerical scheme described in this section has been implemented in CUDA programming language in order to be able to run the model in GPUs.A high efficient GPU based implementation of the numerical scheme allows us to compute high accuracy simulations in very reasonable computational time.
The numerical simulation presented in this work covers a wall clock time period of 10 minutes.14,516 time iterations were necessary to evolve from initial conditions to final state, 10 minutes later.This required a computational time of 1,528.83s (approx 25.5 min), which means 44 millions of computational cells processed per second in a nVidia GTX480 graphic card.

Parameter sensitivity analysis
To overcome the uncertainty inherent to 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).Previously 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 proxy for model sensitivity and will be that m f = 0.001 it is not an admisible 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 Figure 8  In Figure 9, each panel shows the values of the maximum runup computed in each region for a fixed value of the ration 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 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 parameters must be ruled out.Now, in order to find the optimal parameters, according to four criteria that will be defined soon, a more detailed ("microscopic" grid of parameters) is used for performing additional numerical experiments.Taking into account the results presented in Figures 7, 8, and 9, the ratio of densities, r, now will 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, 0.4, 0.42, 0.44, 0.46, 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: C1 the runup on the spur southwest of Gilbert Inlet had to be the closest to the optimal 524 m,

5
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 the Cenotaph Island, opening a narrow channel through the trees Miller (1960), C4 the trimline maximum distance of 1,100 m from 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: Setting the three parameters to the values given above, simulation satisfied the previous four criteria with very good accuracy, SC3 the wave produced a narrow channel crossing through Cenotaph Island (Fig. 10.c), and, finally, SC4 the runup reached more than 1,100 m distance from high-tide shoreline in Fish Lake area (Fig. 10.c).
Figure 10 graphically shows the four selected criteria enumerated above.In Figure 7 the maximum runup values for the optimal solution in the four regions selected for the sensitivity analysis are marked in the four panels.Next section describes in some detail the numerical experiment performed with the optimal set of parameters.

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 is presented: first the main characteristics of the giant wave generated in Gibert 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 s -39 s Following landslide trigger, the generated wave reaches its maximum amplitude, 272.4 m, at t = 8 s (Figure 11.b).The wave spreads outwards in a southward direction with decreasing amplitude (Figures 11.c,11.d).22 s after triggering, the wave hits the bottom of southwest spur of Gilbert Inlet.At t = 39 s the maximum 523.9 m runup is reached (Figure11.f).t: 5 min -8 min Large inundated areas were formed both around Fish Lake, in the north shoreline, as in the flat areas surrounding the Paps, while the main wave reaches the narrow area near La Chausse Spit.During this time, the wave amplitude is larger than 15 m, indeed over 20 m on the north shoreline.At 5 min 50 s the wave reaches La Chausse Spit, passing over before reaching the sea, and then partially reflecting a wave back into the bay (Figure 13.g, 13.h).

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 (Figure 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.

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. 10.a).There is large extent 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.

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 runup in the north shoreline of La Chausse Spit is also found.

Fish Lake area
The agreement between model and observation around the flat areas surrounding Fish Lake is good (Fig. 10.c).Here, the inundation extent includes vicinity areas of Fish Lake under 40 m 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.

South Shoreline
The computed runup underestimates the trimline over steeper slopes in the eastern third of the south shores (Fig 14).Moreover, the numerical model provides mean runup heights 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 Chausse Spit (Fig. 10.d).

East area of The Paps
Around the central third of the south shoreline, the agreement between model and observed trimline is less coherent (Fig. 10.d).
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.

Cenotaph Island
One of the items to be checked in the sensitivity analysis presented in Section 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 Figure 10.c.

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

10.d).
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 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 2 , 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 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's 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 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): ↵ = 16 (13 ), r= 0.3 (0.44), m f = 0.08 (0.08).
The numerical results for this case, that must be compared with the optimal case selected in Fig. 10, are presented in Figure 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.Even, 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.In the other hand, 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.
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.D. Arcas was supported by the National Oceanic and Atmospheric Administration with this report being PMEL contribution No. 4665.

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

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

Figure 5 .S
Figure 5. Sketch of the two-layer model.Relation among h1, h2 and H.

Figure 6 .
Figure 6.Location and spatial extension of the 4 regions considered to compute maximum runup for the sensitivity analysis.
, 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 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 closer to observations numerical results.

Figure 8 .
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 runup at the four regions defined in Fig. 6 for different values the Coulomb friction angle, ↵, (in different colors, see legend) is presented.The vertical bars show the runup variation for values of r ranging from 0.3 to 0.6.

Figure 9 .
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 runup at the four regions defined in Fig. 6 for different values the friction coefficient, m f , (in different colors, see legend) is presented.The vertical bars show the runup variation for values of ↵ ranging from 8 to 16 .

Figure 10 .
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; and (D) Entrance of the bay up to Cenotaph Island.

Figure 13 .
Figure 13.Wave evolution from t = 3 min 30 s until t = 7 min.In the right hand panels the track of the inundated areas is kept.Left panels for South view, right panels for plan view.

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

Figure 15 .
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 the Lituya Bay from Cenotaph Island; (C) Fish Lake and Cenotaph Island area; and (D) Entrance of the bay up to Cenotaph Island.