2.1 Geomorphological and tectonic setting of Lituya Bay

Lituya Bay is a fjord in southeastern Alaska, originating from glacier retreat (Fig. 1a) 10 000 years ago at the beginning of the current interglacial period (Pararas-Carayannis, 1999), resulting in its present T-shape (Fig. 1b). U-shaped slopes are the main features of the bay, with recent terminal moraine deposits of former Tertiary glaciation periods (Pararas-Carayannis, 1999).

At the head of the bay, the slopes featured very steep walls, ranging from 670 to 1030 m a.s.l. and to more than 1800 m a.s.l. in the Fairweather fault area, about 3 km from Crillon Inlet. Lituya Bay has a length of 12 km and a width ranging from 1.2 to 3.3 km (Fig. 1b), while the entrance is about 300 m in width (Fritz et al., 2001). The northern and southern channels on the side with respect to Cenotaph Island in the middle of the bay are about 650 and 1300 m in width (Fig. 1b). The shores mostly consist of rocky beaches. At the entrance of the bay, La Chaussee Spit (Fig. 1b) represents the terminal moraine resulting from the Last Glacial Period (Pararas-Carayannis, 1999).

The Queen Charlotte and Fairweather faults are situated on the western coasts of Canada and Alaska north of Lituya Bay. They are part of the fault system along the boundaries of the Pacific and North American plates (Tocher and Miller, 1959). The Gilbert and Crillon inlets represent the geomorphological expression of the Fairweather Fault.

2.2 The 1958 Lituya Bay tsunami event

Fritz et al. (2001) report that from 1980 to 2000 four (and probably more) big waves could be verified in Lituya Bay. This occurrence is most likely due to the “unique geologic and tectonic setting of the bay” (Fritz et al., 2001). Miller (1960) reports that several tsunami events happened in Lituya Bay (1853, 1936, and 1958), devastating the forest and reaching an inland run-up of over 100 m a.s.l.

In contrast with many other bays or fjords, the numerous manifestations of tsunamis in Lituya bay are due to several factors. These are (a) its recent environment history (a fjord formed by glacier retreat), (b) the fragile geological and tectonic configuration (steep slopes consisting of fractured rock slopes in a very active fault area), (c) the presence of a great amount of water in the bay and a deep seafloor, and (d) its climate conditions, including intense rain events and periodic freezing and thawing (Miller, 1960).

The earthquake on 9 July 1958, featuring a 7.9–8.3 M_w magnitude, occurred along the Fairweather fault. A rockslide collapsed into the bay at Gilbert Inlet (Fig. 1) (Fritz et al., 2001). Horizontal movements of 6.4 m and vertical movements of 1.0 m were estimated here, as reported in Tocher and Miller (1959). Fishermen that experienced the event spoke of about 1 to 4 min of shaking. The earthquake may have triggered the rockslide (Fritz et al., 2009). The impact of the rockslide generated a huge impulse wave, whose maximum run-up (524 m a.s.l.; Figs. 1c, 2b) is the highest ever recorded in history (Fritz et al., 2001). The wave propagation along the bay resulted in forest destruction and ground erosion (Fig. 1b). Miller (1960) hypothesises that the rockslide is a source of the tsunami by evaluating photographs from the slopes at the Gilbert Inlet before and after the event. He estimates the volume of the main rockslide (30×10⁶ m³) and defines the upper scar limit as being about 915 m a.s.l. (Figs. 1c, 2a).

To be able to distinguish this mass movement from gradual processes and ordinary landslides, Pararas Carayannis (1999) classifies it as “subaerial rockfall”, while Miller (1960) describes it as a mixture between a landslide and a rockfall process according to Sharpe (1938) and Varnes (1958).

Before the catastrophic event, two gravel deltas were located in front of Gilbert glacier (Fig. 1b, c). The rockslide propagated with very high velocity (Ward and Day, 2010), hitting a part of the glacier and the gravel deltas (Miller, 1960). After the event, the glacier front showed a vertical wall (Fig. 2a), since 400 m of ice had disintegrated. At the same time, the deltas had been washed away. Additionally, Ward and Day (2010) describe that the water ran upslope as a surge or splash. The second maximum run-up of 208 m a.s.l. has been identified near Mudslide Creek on the southeastern side of Gilbert Inlet (Fig. 1b). The wave reached an inland distance of 1400 m on the plain in front of Fish Lake (Fig. 1b) on the northwestern side of the bay (Ward and Day, 2010).

Two fishermen eyewitnessed a violent disturbance at the mouth of Gilbert Inlet shortly after the first shaking. They confirmed the rockslide as the trigger of the impulse wave (Ward and Day, 2010). The fishermen estimated the wave's amplitude to be about 15–30 m as it impacted Cenotaph Island (Fig. 1b). On their boats, they additionally experienced a short period with wave heights up to a few metres shortly after the initial wave (Ward and Day, 2010).

2.3 Existing studies on the 1958 Lituya Bay tsunami event

As one of the most studied cases for landslide-generated impulse waves, the Lituya Bay tsunami of 1958 has been and still is of great interest to the scientific community concerning the assessment of natural hazards. Different approaches have been used to reproduce the tsunami, covering analytical studies, physical scale tests, and numerical modelling.

A three-dimensional model of Lituya Bay at a 1:1000 scale was created by Wiegel (1964). Thereby, a run-up of about 3 times the water depth was observed on the opposite slope of the sliding source. Fritz et al. (2001) recreated the 1958 Lituya Bay tsunami in the impact area. They simulated the wave formation and run-up in a two-dimensional physical-scale model as a vertical section of Gilbert Inlet at a 1:675 scale. The results confirm the hypothesis of the rockslide as the source of the impulse wave and the observed maximum run-up. A three-dimensional pneumatic generator for landslide-generated impulse waves was applied by Fritz et al. (2009). The 1958 Lituya Bay rockslide was thereby recreated in a three-dimensional model at a 1:400 scale.

Several studies are based on the use of analytical equations, including those derived from physical-scale test experiments. The main goal of these analyses is to reproduce the wave amplitude and the maximum observed run-up. Concerning the maximum impulse wave height, the results obtained are within the range of 94–162 m (Kamphuis and Bowering, 1970; Noda, 1970; Slingerland and Voight, 1979; Huber and Hager, 1997; Fritz et al., 2004, 2009). Hall Jr. and Watts (1953) and Synolakis (1987) matched the maximum run-up of 526 and 493 m a.s.l., assuming a measured impacting wave height of 162 m and a water depth of −122 m as input (Fritz et al., 2001). Slingerland and Voight (1979) confirmed that a wave height of about 160 m is required to recreate the maximum observed wave run-up. Heller and Hager (2010) applied the impulse product parameter to estimate the main impulse wave characteristics in Lituya Bay. Considering a slide impact speed of 92 m s⁻¹ (Körner, 1976), they predicted a wave height of 179 m based on the wave channel geometry.

Starting from the experiments of Fritz et al. (2001), Mader and Gittings (2002), Quecedo et al. (2004), Weiss and Wuennemann (2007), and Basu et al. (2010) also numerically simulated the 1958 Lituya Bay tsunami event, with a focus on the impact area, by applying the Navier–Stokes hydrodynamics code in two dimensions. Mader (2001) applied the Simulating WAves Nearshore (SWAN) code to model distinct wave trigger mechanisms. These studies stated that a landslide led a tsunami wave to flood Lituya Bay in July 1958 but that it had a mobilised water volume 10 times lower than the volume of the collapsed rockslide. Mader and Gittings (2002) simulated the Lituya Bay tsunami with the Navier–Stokes Adaptive Mesh Refinement (AMR) Eulerian compressible hydrodynamic code (SAGE), which resulted in a maximum wave height of 250 m and run-up of 580 m a.s.l. at the opposite slope.

Pastor et al. (2008) applied a coupling model for displacement and pore pressure together with a generalised plasticity model that describes soil behaviour. Wave propagation is evaluated using a depth-integrated model with fluidised soil rheology. The interaction of the slide and the seawater is simulated with a level-set algorithm that tracks the interfaces between air, water, and solid. They computed a slide impact speed of 110 m s⁻¹, leading to a maximum wave amplitude of 226 m. To simulate the tsunami run-up, Weiss et al. (2009) used a hybrid model approach for the movement of deformable bodies in a U-shaped valley (comparable to Gilbert Inlet). They obtained a maximum wave amplitude of 152 m and a maximum run-up of 518 m a.s.l. Basu et al. (2010) applied the drift–flux model implemented in the CFD software Flow-3D to simulate the landslide-generated impulse wave formation in the impact area of Lituya Bay. Assuming an initial void fraction of 40 % for the rockslide material, they predicted a maximum amplitude of 200 m and a maximum run-up height of around 673 m a.s.l.

The two-dimensional representation of Lituya Bay of Fritz et al. (2001) was also used by Schwaiger and Higman (2007), Chunqi et al. (2016), and Xenakis et al. (2017) in the context of an SPH (smoothed particle hydrodynamics) modelling approach. SPH allows for a better representation and simulation of the landslide material collapse process and its impact on the waterbody.

Ward and Day (2010) developed a new “tsunami ball approach” to simulate the impulse wave formation and propagation along the whole Lituya Bay. They predicted a wave height up to 150 m in the impact area and a run-up height of 500 m a.s.l. They considered a dual source for the tsunami event: the subaerial rockslide and a huge amount of subglacial sediments released in the bay after the rockslide impact in the waterbody. The resulting trimline is overestimated by the dual-source approach, but only the subaerial rockslide as an impulse wave trigger was not enough to explain the whole flooded area along the bay. González-Vida et al. (2019) applied a finite-volume Savage–Hutter shallow water coupled numerical model (HySEA). The resulting numerical simulations succeeded in reproducing most of the features of the tsunami event.

3.1 Pre-event bathymetry and topography

Digital data and cartographic material concerning the bathymetry and topography of Lituya Bay dated before and after the tsunami event are available. None of these data are close enough to describe the exact configuration of the bay shortly before 9 July 1958.

The 1926 and 1940 bathymetry surveys (U.S. Coast and Geodesic Survey, 1926, 1942) show that the northeastern limit of Lituya Bay has a U-shaped valley with steep slopes and a wide flat seafloor, which constantly increases its depth until a maximum of −220 m a.s.l. on the southern side with respect to Cenotaph Island (Pararas-Carayannis, 1999) and then decreases in the direction of the sea. In the area close to the bay entrance, the bay floor is −10 m a.s.l. on average.

Miller (1960) was the first after 1958 to describe the bay before the tsunami event (Fig. 1c). He describes the area between Cenotaph Island and Gilbert Inlet as a wide expanse with depths between −150 and −220 m a.s.l. He highlights the presence of two deltas on both sides in front of Gilbert Inlet. Miller (1960) mapped the pre-event topographic and bathymetric contours. In the post-1958 surveys, these areas and deltas are not present (Ward and Day, 2010). The 1969 chart obtained from the 1959 survey (U.S. Coast and Geodetic Survey, 1969) shows a flat bay floor (green zones in Fig. 4a). A ridge divided the bay floor into two sub-basins: a smaller one southeast of Gilbert Head (−156 m maximum depth) and a larger one on the south front of Cenotaph Island (−150 m maximum depth). A third survey published in 1990 (U.S. Coast and Geodetic Survey, 1990) creates the possibility of estimating the sedimentation rate over time. The two charts, the first from 1942 to 1969 and the second from 1969 to 1990, differ completely (see Fig. 8 in Ward and Day, 2010). The northeastern sub-basin in front of Gilbert glacier is filled by sediment. The bay floor decreases constantly between Gilbert Inlet and the basin in front of Cenotaph Island. Ward and Day (2010) estimated that 3×10⁸ m³ of material was discharged into the bay after the tsunami event, filling it until the −130 m depth contour, resulting in a 70 m thick deposit. From the previous considerations, they propose a hypothesis to justify the whole infill of the bay between 1926 and the 1958 tsunami event. Given that (i) the sedimentation rate is assumed constant during the last century; (ii) in 1936 a landslide collapsed into the bay (where the generated wave was a tenth of the size of the 1958 tsunami; Miller, 1960); (iii) the 1958 rockslide contained 3–6×10⁷ m³ of material (10 %–20 % of the total infill volume); and (iv) soil, subsoil, and bedrock have been eroded by the wave (about 4×10⁶ m³; Miller, 1960), they suggested that the remobilised sediment located under Gilbert glacier contributed to infilling the bay floor during the tsunami event. The volume generated from the displacement of the deltas in front of Gilbert glacier has thereby not been considered.

All of these considerations are useful to give a good interpretation of the bay pre-event configuration (Fig. 3b). The bathymetric data used for this study are provided by the National Ocean Service: Hydrographic Surveys with Digital Sounding. In particular, data from Survey ID: H08492 (1959) are used as reference bathymetry after the event since this survey is the closest to the 9 July 1958 tsunami event, as well as data from Survey ID: H04608 (1926) the closest survey prior to the event. The survey from 1926 does not feature a sufficient spatial resolution to provide an acceptable bathymetry in the whole bay. Nevertheless, it provides sufficient information on the pre-tsunami bathymetry in the areas at Gilbert Inlet and south of Cenotaph Island. As mentioned in Ward and Day (2010), the infilling material after the 1958 event covers an area between the −120 and the −130 m depth contours. In the map, contour lines defined by Miller (1960) (dashed red lines in Fig. 3a) from a depth of −122 to −220 m show the bay bathymetry before the infill. From these considerations, the elevation points from the Survey ID: H08492 (1959) are taken as representative of the shallower part of the bay floor (from the surface to −120 m depth). Below −120 m depth, the elevation points from the Survey ID: H04608 (1926) are considered to be representative of the deeper area of the bay floor. Due to limited data, the contour lines defined by Miller (1960) are used to better reproduce the shape of the bay floor. For the flatter areas, lines are set between different elevations points to allow a more accurate interpolation for the bay floor surface reconstruction. The delta in front of Gilbert glacier (Fig. 3a) is reproduced considering the information given by Miller (1960).

The pre-event topographic surface is reproduced starting from the current digital terrain model (5 m resolution) available from the DGGS Elevation Portal of Alaska. Contour lines of 5 m are used to recreate the topography. Contour lines of 1 m are also used to highlight some details that influence the estimation of the flooded area (such as steeper slopes, hills, or specific topographic changes). The observed trimline and the run-up records (red spots in Fig. 3b) are used as references to define the required spatial extent of the model topography.

Manual modifications of the DTM are provided at the Gilbert and Crillon Inlet locations in order to recreate those geomorphological elements that have been displaced and washed away during the landslide impact and wave formation. Crillon Inlet, Gilbert glacier, and the deltas in front of them are recreated starting from the descriptions provided by Miller (1960) and the available cartography (see Fig. 16 in Miller, 1960). Miller (1960) reports a scarred area, pre-dating the tsunami event, located north of the maximum observed run-up (Fig. 2b). He only identifies a little scar exactly under the maximum run-up that could have been formed after the tsunami event. The pre-event shoreline is also reproduced starting from the descriptions provided by Miller (1960).

The bathymetric and topographic surfaces are recreated and exported to stereolithography files (STL) via Rhinoceros 6 software (https://www.rhino3d.com/, last access: 24 March 2020).

3.2 Model setup and computational details

3.2.2 Denser fluid setup

The cliff material consists mostly of amphibole and biotite schists with an estimated density of the undisturbed rock of 2700 kg m⁻³ (Table 1). The sliding mass dimension before the collapse is well known. The thickness of the slide has been defined by Miller (1960). The mass of the rockslide is described as a rock prism with a triangular shape (along a vertical section) with a width varying from 730 to 915 m (Miller, 1960; Slingerland and Voight, 1979; Fig. 1c). The length results in 970 m along the slope (Slingerland and Voight, 1979, Table 1). The maximum thickness results in 92 m; the centre of gravity is located at 610 m elevation (Miller, 1960; Table 1, Fig. 1c). Miller estimated the volume of the sliding mass to be about 30.6×10⁶ m³ with an elevation from 230 to 915 m a.s.l.

Table 1Summary of the governing parameters of the 1958 Lituya Bay tsunami event and related references.

Download Print Version | Download XLSX

Due to the lack of data available, the rheology of the rockslide is difficult to describe. Limited information has recently been provided in the literature (Quecedo et al., 2004; Schwaiger and Higman, 2007; Pastor et al., 2008; Mao et al., 2017; Xenaksis et al., 2017), where artificial rheological parameters are obtained empirically following computational experiments or resulting from back analysis of similar rockslides. Xenaksis et al. (2017) describe the slide phase as a mixture of different materials with different properties (rock, soil, ice, and vegetation), where non-Newtonian shear-thinning rheological properties are assumed. Mao et al. (2017) adopted a rheological model as a non-Newtonian viscoplastic fluid model in order to recreate a deformable fluid-like slide body.

Since in this case there is a lack of directly observed data, and since a proper setup and validation of the rockslide model with the software Flow-3D is possible, a simple concept of a denser fluid is adopted. The bulk slide volume and bulk slide density, 51.0×10⁶ m³ and 1620 kg m⁻³, respectively, are used for the denser fluid simulation (Table 1). As was done in Fritz et al. (2001), the reduced bulk density of 1620 kg m⁻³ considers a void content of 40 % (Table 1). The porosity parameter used is based on data from debris flows observed in the Alpine region (Tognacca, 1999). This is not entirely representative of the real rockslide material but gives an appropriate approximation for rockslide properties.

3.2.3 Model concept analysis

A 3D model of the Lituya Bay topography, idealised as a bucket shape, is assumed for the model concept analysis (Fig. 4a), starting from the information provided by the 2D numerical simulations proposed by Basu et al. (2010), which summarise the experiment of Fritz et al. (2001). The simulation time is set to 60 s. The terrain model and the computational domain are presented in Fig. 4a.

Different slope angles of 35, 40, and 45 ^∘ are set to verify the influence of the impact angle on the impact speed (Fig. 4b). Despite the change in inclination, the difference in heights is maintained at equal levels in all simulations, thus the geometrical setup is done in such a way that this difference is measured from the sea level to the upper edge (915 m a.s.l.) of the denser fluid at its initial position (Fig. 4b). To ensure the same volume and shape, the centre of gravity (as well as the lower edge of the denser fluid) is at different heights for different slope situations.

The computation domain has its origin located on the bay floor, assumed to be −122 m below sea level (Table 2). An orthogonal grid comprising a uniform cell size of 10 m is defined for these models. The grid includes the air space above the bay between the headlands to accommodate the waves according to the VOF algorithm (Basu et al., 2010). In order to reduce the number of active cells in the domain (thus saving computational memory and improving model stability), a solid body occupying higher air cells (where no fluid is expected to occupy space during the simulation) is set as a domain remover. This is applied the same way for the model concept analysis as for the formation and propagation models, in a way that the limits of the active domain correspond to the limits of the recreated topography (Fig. 5).

Table 2Summary of the simulation setup and modelling tasks.

Download Print Version | Download XLSX

The boundaries are specified as the outflow on the free sides of the idealised topography to allow the fluid to flow out of the model without any kind of interaction or reflection. The extent of the flow domain is set in such a way that the fluid interacts mostly with the boundary that represents the bay floor and inland slopes (left and right boundaries).

The presence of the glacier, and possible virtual walls, to constrain the fluid during its movement along the slope, is also considered (Fig. 4a). With regards to the evaluation of the modelling concept, it is expected that the impact velocities stay within the interval 90–110 m s⁻¹ (Table 2).

The initial fluid in the bay represents typical seawater conditions and features a density of 1035 kg m⁻³ (Table 1).

3.2.4 Modelling of the impact area and the whole bay

Further simulations in the impact area (near field) include the topography surface and the recreated bathymetry (Fig. 5, Table 2). The simulation time is 70 s. Different uniform cell sizes are set up for these simulations (20, 10, and 5 m) in order to verify the accuracy of the results based on the grid resolution (Table 2). The simulation domain extends to 1600 m × 4000 m in the x–y direction and 1200 m in elevation. The same boundary conditions as those used for the model concept analysis are set for these simulations. The denser fluid shape is redefined starting from satellite images and cartographic material pre-event. The resulting volume is readapted to the detachment area. The maximum used slide thickness of 134 m is equivalent to 1.4 times the thickness of 92 m provided by Miller (1960). This increase of 40 % in thickness was also considered by Fritz et al. (2001). They adopted this rise in slide density to compensate for the void fraction current in granular flow to match the slide mass-flux per unit width.

For the impulse wave propagation along the whole bay (far field), the domain extends to 6810 m × 13 575 m in the x–y direction and 1200 m in elevation (Fig. 5, Table 2). The simulation takes 7 min. Uniform and non-uniform cells of different sizes are set up 20 m × 20 m × 20 m, 20 m × 20 m × 10 m, and 15 m × 15 m × 10 m). At the domain limits, at the Gilbert and Crillon inlets, and at the seaside, the outflow boundary condition is set to allow the wave to flow out from the model domain.

Control points (Fig. 5) representing specific records of run-up are set in order to validate the results. Several observation gauges (history probes) are set along the entire model domain to achieve information regarding impact time; impact speed (defined as the one of the centre of gravity of the sliding mass entering into the waterbody; Evers et al., 2019); wave propagation speed; and characteristics such as water surface elevation (or wave amplitude), flow speed (magnitude of the velocity vector, resulting from the vector components x, y, and z of the fluid at a specific point in the 3D domain), and their trend in time. In the impact area, probes P1, P2, and P3 are located along the main wave flow direction (Fig. 5) for a streamwise distance from x₀ of 45, 688, 1342 m, respectively, from the impact point. Other history probes are set parallel to the bay length (Fig. 5), starting in front of the delta corresponding to Cascade glacier for a streamwise distance from x₀ of 600 m (P4), 3100 m (P5), 5600 m (P6), 6600 m (P7N/P7S, both located laterally with respect to Cenotaph Island), 8100 m (P8), and 10600 m (P9).

The surface roughness in Flow-3D consists of two components. The first results from the pre-processing phase of the considered solid structures (STL files) with the FAVOR method. Depending on the mesh structure and grid resolution, it features divergences from the original solid structure. The computational geometry usually features a rougher surface than the solid structure in cases where the mesh orientation does not fit perfectly with the surface slope. Furthermore, a second additional roughness parameter, defined as equivalent grain roughness (m), can be set for each solid structure to consider, for example, vegetation.

The computations in this study are mainly set up with the equivalent grain roughness equal to 0 m for the topographic surface. In order to verify the influence of the vegetation on the inundation process and the trimline definition, simulations with values of 1 and 2 m of roughness are set up for simulations with a grid resolution of 15 m × 15 m × 10 m.

Figure 6Impact speed distribution vs. time for the impacting fluid considering different slope angles of (a) 45^∘, (b) 40^∘, and (c) 35^∘. The lower and upper limits represent a reasonable speed interval for the centre of gravity of the deforming fluid when entering the waterbody (from Eq. 3.5 in Evers et al., 2019, with a dynamic bed friction angle δ of 0 and 14^∘).

Download

4.1 Evaluation of the denser fluid concept

Several preliminary simulations are accomplished with the numerical model, in order to test the concept of the denser fluid with respect to the seawater density for the impacting fluid. Different configurations are investigated: (a) different slope angles (35, 40, 45^∘), (b) absence or presence of Gilbert glacier (as a vertical wall of 100 m a.s.l.), and (c) use of virtual walls to constrain the denser fluid during its movement along the slope. This is done to observe the reaction of the wave depending on the changes in these options, in comparison to the simple bucket shape.

The whole process reflects what resulted from the experiment of Fritz et al. (2001), which describes the high velocity of the slide impact process with the following two main steps: (a) the impact of the slide, with the emergence of cavity effects and generation of the impulse wave, and (b) the collapse of the air cavity together with the mixing process of the different components. The formation of a large air cavity after the initial impact is well observed in the computational model.

The results of the model concept analysis are discussed as follows.

• The denser fluid reaches the waterbody within 10–14 s. The maximum impact speed varies between 92–114 m s⁻¹ depending on the modelled impact angle (Fig. 6c). The denser fluid does not act as a non-deformable body during the moving process along the slope. Thus, this deformation has a considerable influence on the impact speed. In Fig. 6, the fluid speed during the impact process is shown for every considered slope angle. An upper and a lower limit (dashed black lines) define a reasonable speed interval for the centre of gravity entering into the waterbody (values obtained from equation 3.5 of Evers et al., 2019, with a dynamic bed friction angle δ of 0 and 14^∘).

• The presence of the virtual constraining walls does not significantly influence the impact speed, although it avoids the mass spreading along the slope during its movement.

• The impulse wave is formed and reaches its maximum height after 9–10 s following the impact, with a wave amplitude ranging between 203–220 m.

• The presence of the constraining walls increases the wave amplitude by 5 to 15 m. This means that the shape of the denser fluid (fluid thickness), as it is constrained by the walls at the impact, influences the wave characteristics more than the impact speed. The presence of the glacier does not influence the wave formation.

• The additional presence of the glacier, together with the constraining walls, affects and increases the impulse wave just before the impact on the opposite headland (20 s after the slide impact). Here, the wave amplitude ranges between 156 and 217 m. It is observed that the wave has no possibility to complete its breaking process, hitting the opposite headland very violently and starting its run-up process along the slope.

• Different maximum run-up values result from different model configurations. They overestimate the observed one, ranging from 570 to 790 m a.s.l. between 36–38 s after the slide impact and 16–17 s after the wave hits the opposite headland. Once the maximum run-up is reached, a backflow of the wave is observed.

• Closer run-up values of 524 m a.s.l. are found for calculations considering the simple bucket shape of the bay, without the presence of the glacier and walls. Considering these two elements in the model, the maximum observed run-up is highly overestimated.

• It is noticed that the use of different order approaches for density evaluation influences the interaction between the two fluids. With the first-order approach, the mixing process leads to a change in density of the denser fluid from 1620 to about 1250 kg m⁻³. With the second-order approach, the density changes to 1400 kg m⁻³. In both cases, a part of the denser fluid runs up a short distance on the opposite slope. The order of the density evaluation approach does not influence the wave characteristics.

The main aim of several authors was to reproduce the impulse wave formation and reach the observed run-up. Provided that the wave run-up could reproduce this value (and consequently the backflow), it could be assumed to have obtained a reliable result and a good reproduction of the 1958 Lituya Bay tsunami event. However, this is not properly correct when taking into consideration the complete run-up process. Actually, the wave did not stop at 524 m a.s.l. but overtopped the hill crest and continued to flow diagonally along the slope to the other side for a distance of about 1 km before re-impacting the sea (Figs. 1c, 2b). An overestimation of the maximum run-up, in these simplified simulations, makes sense for allowing the further overtopping at the hill crest. More “power” is needed to reproduce the phenomena and what has been observed in the whole bay. For the model concept with a topographic surface, the presence of the glacier (and walls to constrain the denser fluid during the movement along the slope) might be necessary to recreate the impulse wave formation and run-up at the head of Lituya Bay.

4.2 Wave formation and run-up – near field

A topographic and bathymetric surface of the impact area is set up and the shape of the denser fluid is readapted to the detachment area (Figs. 5, 7a). What changes here with respect to the model concept analysis is that (a) the slope angle is not constant but ranges from 45^∘ at higher elevation to 35^∘ at the shoreline, and (b) the volume of the seawater is involved in the numerical model since the deltas where not considered in the simplified simulations (about 1.73×10⁶ m³ of seawater compared to 3.34×10⁶ m³ in the model concept analysis). This can have a significant influence regarding the water volume involved in the wave formation and run-up.

Figure 7Fluid model at (a) 0 s, (b) 8 s, and (c) 12 s impacting the sea, coloured by the depth-averaged flow speed in (m s⁻¹) with a range of 0–100 m s⁻¹. Uniform grid resolution of 5 m.

Download

The main aim of this part of the work is to investigate the wave characteristics after the impact of the denser fluid, in order to simulate not only the maximum run-up but also the overtopping process and the flow path along the slope on the other side with respect to the Gilbert Inlet and to recreate the related trimline.

The detachment area, where the rockslide failed, is confined on the left side from the topographic surface, while on the right side two scar channels are presented (Figs. 2a, 7a). These are related to other smaller rockslides that occurred during the earthquake but were not involved in the impulse wave formation (Miller, 1960). For this reason, a constraining wall (invisible in the images, Fig. 7) is set only on the right side with respect to the rockslide described by Miller (1960). A simulation without the wall is also set up to observe the eventual fluid collapse and impact process.

The results obtained (Fig. 8) vary according to the adopted uniform grid resolution (20, 10, 5 m) (Fig. 5). The best fit with the observed maximum run-up in the impact area can be achieved using a uniform grid resolution of 5 m.

Figure 8Wave formation and propagation in the impact area using the second-order approach for the density evaluation. Observation gauges P1, P2, and P3 are set to verify the water surface elevation and flow speed. Their trends are shown in the graphs for different grid resolutions (R: 5, 10, 20 m). More accurate results are obtained using the grid resolution of 5 m (sky-blue line, R5).

Download

Figure 9Impulse wave run-up on the opposite slope. At the 46 s time step the wave reaches the maximum observed elevation of 524 m a.s.l. (flow depth contours). A part of the wave body overtops the hill and proceeds on its path in a diagonal direction with respect to the slope gradient (shown by the different shades of purple for the time steps of 48, 50, 52, 54, and 56 s).

Download

A description of the wave formation and run-up resulting from the simulation with 5 m grid resolution and the adopted second-order approach for the density evaluation is provided below. This model takes 30 h to run.

• 0–15 s. The denser fluid reaches the sea after 10 s with a maximum impact speed of 93 m s⁻¹ and a maximum thickness of 79 m (P1 – x₀=45 m, Figs. 7c, 8a). The depth-averaged flow speed varies from 40 m s⁻¹ in the upper part to 90 m s⁻¹ in the lower part of the fluid during the movement (Fig. 7).

• 15–30 s. After 24 s after the release of denser fluid, the maximum estimated wave amplitude results in 208 m with a flow speed of 78 m s⁻¹ (P2 – x₀=688 m, Fig. 7b). Slightly further (x₀=885 m), the wave maintains its height to start its breaking process. A part of the wave also flows on the glacier. The wave front runs up the delta and the following slope (Fig. 7c).

• 30–35 s. The whole wave crashes on the opposite headland after 30 s following the denser fluid release (20 s after the impact), with a variable water surface elevation of 129–147 m a.s.l and a flow speed between 50–70 m s⁻¹ (P3 – x₀=1342; Fig. 7c). The wave-breaking stage is not complete, it partially breaks when it flows on the delta.

• 35–50 s. The wave runs up towards the headland and the scars located upon the delta. The maximum observed run-up (524 m a.s.l.) is reached after 46 s (36 s after the denser fluid impact, Fig. 9) with a flow depth of 4–10 m and a flow speed of 12 m s⁻¹. A part of the wave body overtops the hill crest, but a backflow is also observed.

• 50–70 s. The wave flows in a diagonal direction relative to the slope, with a depth-averaged flow speed of 50–70 m s⁻¹. The wave reaches the seaside 8 s after the maximum run-up (54 s from the denser fluid release, Fig. 9). The flow depth is about 15 m with a flow speed of 70 m s⁻¹.

During the collapse process, it is noticed that the left part of the denser fluid is well constrained by the actual topography (Fig. 7b, c). Avoiding the wall on the right side, the mass largely spreads and collapses on the glacier, losing a great amount of volume involved in the impact process and decreasing the wave formation. The presence of the wall constrains the denser fluid on this side and allows it to collapse in the waterbody. In addition, Gilbert glacier acts as a constraining wall, and the delta in front of the glacier acts as a ramp.

The maximum wave amplitude of 208 m is located upon the terminal front of the delta on the bay floor (where the history probe P2 is located; graph in Fig. 8b). Here the wave celerity is estimated to be 55 m s⁻¹ (from Eq. 2.2; provided by Heller et al., 2009). The wave starts to break because of its interaction with the decreasing bay floor depth. Fritz et al. (2001) observed the maximum wave height (>200 m) at x₀=600, while at x₀=885 they reconstructed a wave height of 152 m.

The presence of the scar area on the right side of the maximum run-up has a key role in the run-up process, since it allows the wave to run up along a channel (Fig. 9) and reach the elevation of 524 m a.s.l. This observation supports the topography description provided by Miller (1960). Despite the reproduction of the expected overtopping over the hill and the flow process on the other side of the slope, the resulting trimline appears to be underestimated compared to the observed one (light blue in Fig. 12b).

Additionally, if 524 m a.s.l. is the maximum run-up elevation observed from the trimline, it is not clear in the scar area on the right side, since there is no evidence of a forest trimline. With regards to the simulation results, it appears that the maximum run-up reached an elevation up to 600 m a.s.l. in this part of the slope.

4.3 Impulse wave propagation – far field

The aim of these simulations is to reproduce the wave propagation along the bay, to understand how the waves interact and inundate the inland area, and to recreate the actual trimline. For the wave propagation, only the second-order approach for the density evaluation has been used. Observation gauges for water level measurement allow more insights into the wave characteristics during the propagation and to observe the wave attenuation along the bay to the seaside (Fig. 5).

The results obtained from the simulations vary depending on the resolution of the computational grid. A description of the wave propagation and inundation resulting from the simulation with 15 m × 15 m × 10 m grid resolution is provided (details are referred to the primary wave front). The model takes 165 h to run completely.

• 0–60 s. Over the impact area, the wave starts to propagate, resulting in a water surface elevation of 40 m a.s.l. and a flow speed of 9 m s⁻¹ at P4 (x₀=600, Fig. 10a).

• 60–120 s. The wave propagates in open water; at P5 (x₀=3100, Fig. 10b) a water surface elevation of 39 m a.s.l. and a flow speed of 19 m s⁻¹ are recorded due to the amount of water flowing down the slope with high velocity. An impact on the southern side of the bay in front of the Gilbert Inlet is observed, where a secondary wave front is generated due to reflection. The primary wave front reaches the Mudslide Creek delta and floods the inland area with a depth-averaged flow speed of 20–30 m s⁻¹. The second highest run-up results are overestimated with 233 m a.s.l. about 94 s after the release of the denser fluid (the observed one is 208 m a.s.l.).

• 120–200 s. The wave splits into two fronts approaching and impacting Cenotaph Island (water surface elevation of 22 m a.s.l. and 8 m s⁻¹ flow speed recorded at P6, x₀=5600, Fig. 10c). On the southern side, the wave slightly slows down due to the attenuation process, with a water surface elevation of 19 m a.s.l. and a flow speed of 5 m s⁻¹ at P7S (x₀=6600). The steep slopes on the southern side of the bay are completely flooded by the wave (Fig. 10d). On the northern side of the island, where the bay floor gets more shallow (depth 20–40 m) and narrow, the water surface elevation results in 15 m a.s.l. with a speed of 7 m s⁻¹ at P7N (x₀=6600, Fig. 12d). This is due to a breaking process.

• 200–280 s. Due to diffraction, the waves turn around the island and flood the western side of Cenotaph Island. The wave front from the southern channel comes first, as observed at P8 (x₀=8100, Fig. 10d and e) resulting in a flow speed of 4 m s⁻¹ and a water surface elevation of 10 m a.s.l. The flatter northern side of the bay is flooded (Fig. 10e).

• 280–340 s. The wave reaches the maximum distance of 1400 m, flooding the area in front of Fish Lake with a depth-averaged flow speed of 10–25 m s⁻¹ and a corresponding water surface of 15–5 m from the ground (Fig. 10f). The wave approaches the mouth of the bay, resulting in a water surface elevation of 13 m a.s.l and a flow speed of 5 m s⁻¹ at P9 (x₀=10 600, Fig. 10f). The second wave front reaches the first one, resulting in a long period wave; it takes 180 s to pass over the history probe P9 (from 240 s to the time limit of the simulation, 420 s).

• 340–420 s. After 340 s following the release of the denser fluid, the wave reaches the seaside, completely flooding La Chaussee Spit and the nearby areas with a depth-averaged flow speed of 10–20 m s⁻¹.

Figure 10After the fluid impacts into the sea, the wave propagates and floods the inner land along the bay. The images show the flow speed at (a) 40 s, (b) 70 s, (c) 140 s, (d) 200 s, (e) 240 s, and (f) 340 s. Different observation gauges are set to check the wave attenuation during wave propagation. The trend of the water surface elevation and flow speed are shown in the graphs for the related observation gauge, adopting different non-uniform grid resolutions (R: 20 m × 20 m × 20 m, 20 m × 20 m × 10 m, and 15 m × 15 m × 10 m). The purple colour of the inland are represents the inundated areas when the wave propagates (the simulation video is listed in the Data Availability section below).

Download

The main wave is mainly responsible for the forest destruction, but secondary reflected waves along the bay also contribute to the observed trimline. A clear example is the wave reflected from Mudslide Creek impacting the opposite northern slope of the bay at 140 s (Fig. 10c). Other secondary wave fronts spread from the bay head due to several reflections of the backflow in front of the Gilbert Inlet.

In Fig. 11a, the wave propagation for the primary wave front is illustrated. The reported values represent a mean value of wave propagation speed (from 40 to 17 m s⁻¹) for each space interval, starting from the records provided by the gauges located along the bay, considering the wave front position at the time of its passage upon every singular gauge (Fig. 11b) and adopting flow path lines for distance estimation. These reflect the interaction of the impulse wave with the bay floor (thus local changes in water depth) and the bay shape. Higher values of the mean propagation speed are estimated between x₀ (impact location) and P8, ranging between 36–40 m s⁻¹, and vary in response to the local change of the water depth (−150–220 m), inducing slight local increases or decreases in the mean propagation speed. A value of about 60 m s⁻¹ is observed between P4 and P5, probably due to the water flowing down the slope, whose impact into the sea induces a local acceleration of the wave in open water, and thus a higher propagation speed. Cenotaph Island and the northern shallower bay floor represent obstacles that induce a slowdown for the wave propagation speed. Due to the sudden shallowing of the bay floor (about −20 m in water depth) wave speed decreases drastically between P6 and P7N (28 m s⁻¹). A general decrease in the propagation speed is observed from P7S to the mouth of the bay (from 33 to 17 m s⁻¹), related to the progressive decrease in the water depth (from –160 to –10 m). The wave attenuation process, in terms of wave amplitude (from 40 to 18 m), proceeds from the head of the bay until the seaside. Dashed lines represent the secondary wave in time. Its role becomes relevant after gauge P6: the second front approaches the first one evolving in a whole wave body between P9 and La Chaussee Spit (C.S. in the graph), inducing an increase in the wave amplitude (from 13 to 18 m) before a breaking process (due to the high shallowing of the bay floor compared to the wave amplitude).

Figure 11(a) The wave attenuation process during its propagation in the bay. The first and the second wave fronts are represented by the full and dashed lines respectively. Mean wave propagation speed is estimated starting from the records at every gauge, considering the position of the first wave front displayed in (b) at the moment of its passage upon the gauges and path lines to estimate the distance of propagation (background topography from © Google Earth Pro 7.3.2.5776; last access: 24 April 2020).

Independently from the grid resolution, general discrepancies in the trimline definition are observed (Fig. 12a, b). Some results are underestimated, as is the case, for example, on the slopes on the southern side of the bay head, the western part of Crillon Inlet, and the Mudslide Creek location. Others are overestimated, such as at the Cascade glacier location, the second highest run-up after the Mudslide Creek, and southern part of Cenotaph Island.

Figure 12(a) Different results of the inundated area and the related trimline, with respect to the observed data (red line), for different grid resolutions and relative roughness equal to 0 m. Sections are reported in Fig. 13. (b) At Gilbert Inlet, the resulting trimlines are defined from the grid resolutions used for the impact area simulations (20, 10, 5 m) and relative roughness is equal to 0 m. The dashed yellow line represents the shoreline before the tsunami event. (c) The resulting inundation area varies depending on the selected relative roughness of 0, 1, or 2 m for the topographic surface (background topography from © Google Earth Pro 7.3.2.5776; last access: 24 April 2020).

The adoption of different values of relative roughness for the topographic surface (0, 1, 2 m, Fig. 12c) results in an evident change for the inundation process. As shown in Fig. 12c, important differences in the flooded area are evident on flatter locations, mainly present in the western region of the bay. Additionally, adopting a roughness of 2 m, the second maximum run-up at the Mudslide Creek results in 210 m a.s.l. Therefore, the trimline obtained from the simulation, with 2 m of relative roughness for a grid resolution of 15 m × 15 m × 10 m, is very close to the one observed, even if some small underestimations and overestimations are still present. However, it is noticed that the simulation with a uniform grid resolution of 20 m can also reproduce the tsunami trimline well.

The fluid mixture process and the submerged propagation of the denser fluid along the bay floor takes place using the second-order approach for the density evaluation. At the end of the simulation, the denser fluid reaches a distance up to 4 km from the impact point, still propagating with a low speed of 5 m s⁻¹ and a thickness of about 35 m. The bulk slide density of the denser fluid decreases during the propagation from 1620 kg m⁻³ to approximately 1080 kg m⁻³, which is close to the seawater density.