## Abstract

Tsunamis caused by landslides may result in significant destruction of the surroundings with both societal and industrial impact. The 1958 Lituya Bay landslide and tsunami is a recent and well-documented terrestrial landslide generating a tsunami with a run-up of 524 m. Although recent computational techniques have shown good performance in the estimation of the run-up height, they fail to capture all the physical processes, in particular, the landslide-entry profile and interaction with the water. Smoothed particle hydrodynamics (SPH) is a versatile numerical technique for describing free-surface and multi-phase flows, particularly those that exhibit highly nonlinear deformation in landslide-generated tsunamis. In the current work, the novel multi-phase incompressible SPH method with shifting is applied to the Lituya Bay tsunami and landslide and is the first methodology able to reproduce realistically both the run-up and landslide-entry as documented in a benchmark experiment. The method is the first paper to develop a realistic implementation of the physics that in addition to the non-Newtonian rheology of the landslide includes turbulence in the water phase and soil saturation. Sensitivity to the experimental initial conditions is also considered. This work demonstrates the ability of the proposed method in modelling challenging environmental multi-phase, non-Newtonian and turbulent flows.

## 1. Introduction

The Lituya Bay tsunami and landslide, which occurred in Alaska in 1958, was triggered by an 8.3 Richter magnitude earthquake, resulting in an estimated 3×10^{7} m^{3} of the rock sliding into a bay generating a wave run-up of 524 m and partial overtopping [1]. Figure 1 shows an aerial photograph of the Lituya Bay fjord overlaid with a graphical representation of the landslide phase and the run-up area.

Scientific interest has gathered around this event due to the devastating consequences that similar phenomena may have. A variety of different computational studies have been performed to reconstruct this event (e.g. [4–9]), with notable examples being the work of Basu *et al.* [3], who used commercial volume-of-fluid software, the smoothed particle hydrodynamics (SPH) simulation of Schwaiger & Higman [10] and the Sanchez-Linares *et al.* [11] Monte Carlo finite volume method for shallow water equations. All of these studies have used the experimental work of Fritz *et al.* [2] as a reference. Among the computational methods mentioned above, the ‘weakly compressible’ SPH (WCSPH) of Schwaiger & Higman [10] is found to exhibit the closest agreement with the experiment.

It is generally observed that although the aforementioned methodologies manage to reproduce accurately the wave run-up height reported in the experiment of Fritz *et al.* [2], they either unsuccessfully reproduce the landslide-entry profile or fail to present it altogether. The scope of the current study is to reproduce accurately the physical modelling of the landslide and wave run-up process by applying and generalizing the incompressible SPH method to include appropriate physics.

SPH, first applied in hydrodynamic flows by Monaghan [12] and Takeda *et al.* [13], is a versatile computational method for studying flows involving free surfaces and multi-phase effects, since it is based on a fully Lagrangian discretization of different fluid domains where fluid particles of constant mass are used as interpolation points. SPH has been successfully used in the past in problems which involve landslide-generated waves (e.g. [14–18]), proving the suitability of the technique for such applications. For a comprehensive review of numerical models of landslide-generated tsunamis, see Yavari-Ramshe & Ataie-Ashtiani [19]. In the current work, the incompressible SPH (ISPH) with the shifting methodology of Lind *et al.* [20] and its non-Newtonian extension due to Xenakis *et al.* [21] are used, since both were found to perform well for free-surface incompressible flows in providing a good representation of the evolution of the flow, including pressure distributions, for Newtonian and inelastic non-Newtonian flows.

The work presented herein shows that ISPH gives a good representation of the Lituya Bay tsunami and landslide, given careful consideration of the physical phenomena involved and the initial conditions. The idealized experiments [2] related to the 1958 Lituya Bay tsunami and landslide are reproduced. It is initially shown that, with a preliminary ISPH method which follows the physics used in the computational model of Schwaiger & Higman [10], the experimental wave run-up height can be reproduced but the landslide-entry profile is not well predicted. The ISPH model is then extended to include a turbulence model for the water phase and a saturation model for the landslide phase. The experimental initial conditions, where a pneumatic landslide acceleration was imposed, are also implemented. The final results show a close agreement with the benchmark experimental findings [2], validating a novel and more complete computational methodology applicable to complex geotechnical phenomena as embodied in the Lituya Bay landslide and tsunami.

## 2. Method

### (a) Governing equations

The mass and momentum conservation equations for incompressible flows are as follows:

where ** u** is the velocity field,

*ρ*is the density,

*p*is the pressure,

**is the stress tensor and**

*τ***F**the body force. In this study, both laminar and turbulent regimes are considered, the latter within a Reynolds averaged Navier Stokes (RANS) formulation for turbulence modelling. Any flow variable

*A*is now written in terms of the statistical mean

*A*′ represents the fluctuations of the field variable

*A*[22]. Taking into account the Boussinesq approximation for the Reynolds stress tensor, the RANS equations can be written as:

*k*being the turbulent kinetic energy. For Newtonian flows

**D**being the shear strain tensor (shown for two dimensions in equation (2.18)), while

*μ*

_{E}=

*μ*+

*μ*

_{T}with

*μ*being the laminar viscosity and

*μ*

_{T}the eddy viscosity. In the rest of the text, the RANS notation presented in equation (2.2) will be omitted for simplicity, but note that in turbulent regimes the Reynolds-averaged notation of equation (2.2) is implied.

### (b) Divergence-free incompressible smoothed particle hydrodynamics with shifting

The computational method used in this work is based on the divergence-free ISPH method [23], in which any variable can be approximated as an interpolation over the computational domain:
*V* represents the computational domain and the smoothing length, *h*, represents the effective width of the kernel *W*. In the SPH formalism, equation (2.3) is discretized to give:
*A*_{j}, *m*_{j} and *ρ*_{j} are the value of any field variable *A*, the mass and the density of particle *j* at position **r**_{j}, respectively. The summation is performed over particles which lie within a circle of a radius determined by the kernel *W* centred at **r**_{i}. A quintic spline kernel [24] is used with a smoothing length *h*=1.3 d*x*, where d*x* is the initial particle spacing.

In the divergence-free methodology proposed by Cummins & Rudman [23], the particle positions, *n*, are advected with velocity *t* is the time-step size. At these positions, an intermediate velocity field, *n*+1:
*W*_{ij}=*W*(**r**_{i}−**r**_{j},*h*) and *η* is a small number introduced to avoid a singularity in the denominator. By substituting equation (2.8) into equation (2.7), a linear system is retrieved which is iteratively solved using a bi-conjugate gradient stabilized (BiCGSTAB) solver. Further details can be obtained in Xu [26]. The pressure gradient is next added to *n*+1 time-step according to
*et al.* [20] and Xenakis *et al.* [21], respectively, while it has recently proved to be beneficial in multi-phase SPH simulations (e.g. [29,30]). In this work, the compound shifting method [31] is used with
*δ**r*_{s} is the shifting distance, *A* is a problem-independent dimensionless constant with a value of 2 [28], and

#### (i) Turbulence

The divergence-free ISPH method with shifting is also implemented when turbulence modelling is considered. This method is now used to discretize and solve the governing equations for RANS turbulence models as expressed in (2.2), specifically for the *k*−*ϵ* turbulence model [32], giving the Reynolds-averaged approximation of the turbulent flow field. The *k*−*ϵ* model has been widely applied and validated both in WCSPH (e.g. [33–36]) and in ISPH with shifting [37]. In the current work, the *k*−*ϵ* model is applied along with an improved implementation of the dummy-particle boundary conditions [38] and, as explained in §2g, is validated against direct numerical simulation results for turbulent Pouiselle flow and applied to the more general case of turbulent fish pass flow (presented in detail in [31]).

The evolution in time of the turbulence kinetic energy *k* and the turbulent dissipation *ϵ* reads in the SPH formalism as
*P*_{i} being a production term [33,37], equal to:
*μ*_{T} is related to the turbulent kinetic energy and dissipation as
*μ* and used to calculate the intermediate velocity field **u***. Terms *C*_{μ}, *C*_{ϵ,1}, *C*_{ϵ,2}, *σ*_{k} and *σ*_{ϵ} are constants specific to the *k*−*ϵ* model, taking the standard values presented in table 1. It is noted that equations (2.12) are solved explicitly with the right-hand side computed using values from time step *n*. More details on turbulence modelling in SPH can be found in Violeau [33,35], Leroy *et al.* [37] and Ferrand *et al.* [36].

### (c) Viscous terms

In this work, the divergence of the stress tensor as presented by Shao & Lo [25] is used
** τ** for inelastic non-Newtonian flows is given by the following constitutive law [39]:

**D**| is the second principal invariant of the shear strain rate,

**D**=∇

**+∇**

*u*

*u*^{T}, [39]:

*μ*

_{nN}is the effective laminar viscosity of a non-Newtonian phase, calculated as a function of |

**D**| in §2d. Equation (2.16) is simplified to

**=**

*τ**μ*

**D**for Newtonian flows. In two dimensions, the shear strain rate tensor can be written as

*D*

_{ij}, of tensor

**D**are obtained using finite difference approximations, before decomposing into

*x*- and

*y*-directions [25]. Thus,

*μ*

_{E}=

*μ*=

*μ*

_{nN}.

### (d) Non-Newtonian rheological models

In this work, three inelastic non-Newtonian rheological models [39] are considered.

*Power-law model*: The effective viscosity in a power law model is
*μ*_{pl} and *N* are the fluid consistency coefficient and the flow behaviour index, respectively.

*Bingham model*: The Bingham model can be represented by a multi-valued function, where the fluid does not deform when the shear stress is below the yield stress *τ*_{Y} (thus resembling solid material behaviour). When the shear stress exceeds *τ*_{Y}, the fluid adopts a constant viscosity *μ*_{B}, thereby resembling a Newtonian fluid. Therefore, the constitutive equation of a Bingham fluid effectively becomes:
*et al.* [41] is implemented, with
*α* is a large scalar (order of 10^{2}).

*Herschel–Bulkley model*:
*μ*_{HB} is the fluid consistency coefficient. As with the Bingham model, in the Herschel–Bulkley rheological model the solid zone could be replaced by a highly viscous fluid or a function which approximates the Herschel–Bulkley constitutive equation. In this study, a bilinear model similar to equations (2.22) was used for the Herschel–Bulkley simulations.

### (e) Saturation model

Saturation describes the process through which voids in a porous material are filled with a liquid, typically water [43]. Such phenomena are crucial in many environmental and industrial applications where a sedimentary material (in this case the landslide) becomes saturated due to its interaction with water, influencing its flow characteristics [43,44]. In the current study, a simple saturation model is considered for the landslide with results shown in §3d(iii).

Saturation models have previously been implemented in SPH applications (e.g. [30,45–48]). The seepage force and the water pressure along with yield criteria, such as the Mohr–Coulomb and the Drucker–Prager, are commonly implemented in the aforementioned techniques. By contrast, in the current work a much simpler model is implemented based on the findings of Mitarai & Nakanishi [44], who have described the rheological transition from dry to wet sediment material as a transition from shear-thinning to shear-thickening behaviour. This transition can clearly have a dramatic effect on the behaviour of the material undergoing saturation. A transition period *T* is assumed in this work, during which the landslide material is saturated. This period although constant, is initialized for each SPH landslide particle separately, when the initial water level is reached. The parameters that change over this time period are the yield stress *τ*_{Y}, the viscosity *μ*, the power-law exponent *N* of the Herschel–Bulkley model and the density *ρ*. A linear transition is considered for each of these parameters which read as
*N*_{i} of particle *i*, with *N*^{0} and *i* reaches the initial water level. The full time history of the rheological parameters is given by

### (f) Multi-phase treatment

To allow the implementation of multiple-phase flows, a colour function is introduced:
*i* [49,50], enabling the different fluid phases to be identified. Moreover, the sharp variations of density *ρ* and viscosity *μ* across the interface can be avoided with the use of a smoothed colour function [49,50]
*ρ*_{a} and *μ*_{a} are the physical density and viscosity of phase *a*, which remain constant through the simulation [50,51]. With this procedure, the continuity of stresses is also satisfied since on the interface, which is defined from

Similar to Mokos *et al.* [29], the normal-to-interface shifting is restricted for one of the two phases. This has been found to minimize any artificial numerical mixing of the two phases, caused by the shifting of different phase particles across the interface.

### (g) Boundary conditions

A number of different boundary conditions are used for the application examined in this work, all of which are based on the dummy boundary particle conditions of Koshizuka *et al.* [38]. In addition to the traditional dummy particle boundary conditions for no-slip boundaries [38], appropriate modifications are made to accommodate free-slip boundary conditions (discussed in [31]) and turbulent flows, where the log law is used to determine particle velocity for dimensionless wall distance *y*^{+}>30. The latter improved upon traditional turbulent SPH implementations using dummy-particles [33,34] through the introduction of a piecewise interaction between boundary and fluid particles (discussed in detail in [31]).

## 3. Results and discussion

### (a) Introduction

A detailed analysis of the Lituya Bay tsunami and landslide results is now given. Initially, the computational problem is defined and preliminary computational results are shown which, following the method of Schwaiger & Higman [10], include a laminar approach for both phases and a gravitational acceleration for the landslide which retains its rheological characteristics without considering any saturation effects. This preliminary analysis shows that although ISPH with shifting and non-Newtonian rheology can improve upon the results of Schwaiger & Higman [10] it fails to capture both the wave run up and landslide-entry profile presented in the benchmark experimental work of Fritz *et al.* [2]. This motivates further investigations into the physical phenomena during the landslide-tsunami process, which leads to the introduction of a RANS turbulence model for the water phase and a saturation model for the landslide. Moreover, the initial experimental conditions are also carefully reproduced. The final results show that the proposed ISPH method can reproduce fully the complex geological events of the Lituya Bay landslide and tsunami as presented in the experiment [2], which shows that the turbulent and saturation phenomena are important physical processes in the problem.

### (b) Problem definition

Figure 2 shows the configuration of the computational domain [2,3], which involves a two-dimensional prismatic channel, with 1000 m height at either side on a 45° angle, separated by a water channel of 122 m depth and 1342 m length.

In the current work, two different input conditions are considered for the landslide. In the first one, a circular-segment-shaped mass of landslide is lying on the North East slope with a width of 970 m and height of 92 m. The second approach corresponds to the experimental initial condition, where a pneumatic accelerator mechanism was used to propel the landslide towards the mass of water. Details of this mechanism and its computational modelling are shown in §3d(ii).

The two media simulated in this multi-phase problem are water and a non-Newtonian phase representing the landslide. Typical values of density and viscosity are chosen for the water, with *ρ*_{water}=1000 kg m^{−3} and viscosity *μ*_{water}=0.001 Pa⋅s, while for the landslide phase a density of *ρ*_{landslide}=1650 kg m^{−3} is chosen, as reported in the experiment of Fritz *et al.* [2]. The whole domain is discretized with an initial particle spacing of d*x*=6 m corresponding to a total of 9219 SPH particles. A challenge for this validation case is to determine accurately the rheological parameters of the landslide phase, as such details are not provided in the literature. In reality, the landslide phase would consist of a mixture of materials with different properties (e.g. rocks, vegetation, soil and ice). Moreover, the landslide phase is expected to have non-Newtonian shear-thinning rheological properties. To determine the rheological parameters of the landslide phase, computational experiments were performed with the Bingham, power-law and Herschel–Bulkley models as presented in §2d.

### (c) Initial comparisons with the weakly compressible smoothed particle hydrodynamics method of Schwaiger & Higman [10]

Initially, a laminar approach is considered for both the landslide and the water phases (following [10]). Moreover, the landslide is accelerated only under the influence of gravity.

The Bingham rheological model (equations (2.22)), the power-law (equation (2.20)) and the Herschel–Bulkley model (equations (2.23)) have been tested, for which, as shown in Xenakis [31], the Herschel–Bulkley model shows the best performance, for *μ*=10^{5.45} Pa⋅s, *τ*_{Y}=100 Pa and *N*=0.375 for a landslide-phase density of 1650 kg m^{−3}.

Comparisons are performed against the WCSPH computational results of Schwaiger & Higman [10] who modelled the initial profile of the landslide phase with a wedge-like shape unlike the circular segment profile used in other published work [3,8] and in the current work. Moreover, the landslide phase was assumed to have a Newtonian rheological behaviour, while the water phase was modelled as an inviscid fluid. Figure 3 shows the results of the ISPH method with the parameters mentioned before against the WCSPH results by Schwaiger & Higman [10]. As shown the non-Newtonian rheology with the ISPH offers improved results compared with the WCSPH Newtonian approach. The improvement of the wave height at position *x*=885 m (figure 3*b*) is notable, where ISPH achieves a closer agreement with the experimental results throughout the computational simulation. Additionally, improvements are shown in the prediction of the landslide-entry profile, which is measured at *x*=−67 m (figure 3*c*) and the wave run-up (figure 3*a*). Despite the improved results shown in figure 3*c* by ISPH, there are still considerable differences comparing with the experimental data. This large deviation is because of the different input methods discussed previously, with Fritz *et al.* [2] using a pneumatic accelerator mechanism (presented in detail in [52]) to propel the landslide mass towards the water, while in the WCSPH method [10] and the current method a gravitational acceleration was acting solely on the SPH particles.

It is interesting to note at this point that the initial shape of the landslide appears to play an insignificant role in the amplitude and length of the resulting wave, which is strongly affected by the shape and velocity of the landslide at impact. Other modelling parameters, such as the rheological behaviour of the landslide phase, and its interaction with the boundary wall appear to have a more direct effect on the wave run-up height, as shown in the parametric study of Schwaiger & Higman [10] and also observed in the parametric analysis performed as part of this study. In fact, there are a number of other computational studies which consider different initial shapes for the landslide with small differences in the generated wave and run-up height observed (e.g. [3,7,10]) but generally poor agreement shown with the experimental profile at the entry phase.

In the following paragraphs, new modelling approaches are introduced which reduce the differences observed in the entry phase of the landslide. The experimental initial conditions [2,52] are reconstructed, to reproduce the landslide-entry profile. Modelling of turbulence and saturation for the water and landslide phases, respectively, is also included.

### (d) The final incompressible smoothed particle hydrodynamics approach

#### (i) Introduction

As shown in §3b, the current state-of-the-art approach can reproduce the wave run up, but generally fails to capture the landslide-entry profile, because of the different initial conditions used in the experimental work of Fritz *et al.* [2] and the numerical simulations. The pneumatic accelerator used in [2] is thus reconstructed using SPH, based on the dimensions presented in Fritz & Moser [52]. Additionally, turbulence is considered for the water phase and saturation for the landslide, as discussed in detail in the following paragraphs.

#### (ii) Reconstruction of the experimental initial conditions

In the current section, the pneumatic accelerator presented in the experimental configuration by Fritz & Moser [52] is computationally reconstructed using SPH particles. Before moving to the modelling of the accelerator the scale of the problem has to be considered. As shown in §3b, the Lituya Bay landslide and tsunami were originally modelled in the actual Lituya Bay fjord’s scale. However, since we wish to make a detailed comparison with the experiments of Fritz *et al.* [2] we model on the experimental scale (i.e. 1:675), which also permits the application of the *k*−*ϵ* turbulence model with a smaller number of particles. For consistency with the experiment [2] and the results shown previously (see §3b) the results shown in the current section are presented in the original Lituya Bay scale, using a Froude similarity model [53].

Figure 4 shows the dimensions of the pneumatic accelerator as presented in [52]. This is reconstructed with SPH particles as shown in figure 5. As illustrated, the flow is accelerated with free-slip boundary conditions on the sloping boundary, as an approximation to uncertain low friction rolling conditions. Notably, Basu *et al.* [3] used free-slip boundary conditions on the landslide wall, but did not present comparisons with experimental data. Moreover, it was found that no-slip boundary conditions result in a trailing effect of the landslide phase and an overall poor agreement with the experimental data as shown in Xenakis [31].

The velocity profile of the pneumatic mechanism was approximated with a combination of two sinusoidal functions as
*u*_{Piston} is the velocity of the pneumatic accelerator, *U*_{Piston} the maximum velocity of the pneumatic accelerator and ^{−1} [2] was essential to approximate the landslide-entry profile.

Results of this implementation are presented in the following section (see §3d(iii)) for rheological parameters which have been determined empirically following computational experiments. As shown, the reconstruction of the landslide pneumatic accelerator has proven crucial in the accurate reproduction of the landslide-entry profile.

#### (iii) Final incompressible smoothed particle hydrodynamics results

In the current section, the final ISPH results are presented which include saturation of the landslide phase and turbulence modelling in the water phase, as outlined in §2. Moreover, the computationally reconstructed pneumatic accelerator is implemented, allowing for the accurate modelling of the experimental initial conditions.

Similar to the laminar case of §3b, a parametric analysis is performed in order to determine the non-Newtonian characteristics of the landslide in this case both for dry and wet states. After examining the Bingham (equations (2.22)), the power-law (equation (2.20)) and the Herschel–Bulkley models (equations (2.23)), the latter was chosen with the rheological values presented in table 2 and for a transition period of *T*=1.95 s in equation (2.25), since it matched most closely with the experimental data [2]. The change of the density corresponds to a 39% void-fraction [2], assumed fully filled with water during saturation. Moreover, the transition of the rheological parameters are based on the findings of Mitarai & Nakanishi [44], who have reported that the rheological behaviour of a sediment shows a transition to shear-thickening behaviour during saturation. In the saturation model used herein a transition from a shear-thinning Herschel–Bulkley model to a pseudo-plastic Bingham model, with higher yield stress and viscosity, was found to give the closest results to the experiment for the wave run-up. Figure 7 shows the shear-rate/shear-stress (normalized with yield stress) relation for the two different rheological models corresponding to the dry and saturated sediment conditions, showing clearly the shear-thinning and pseudo-plastic rheological behaviours.

It should be noted that with the rheological properties shown in table 2 the Reynolds number of the landslide phase is found to be *Re*≈60, which is well within the laminar range. Therefore, a laminar model as defined by equations (2.1) is implemented. The *k*−*ϵ* model is then applied solely to the water-phase along with the saturation model in the landslide phase, since implementation of turbulence alone was found to improve marginally the overshooting wave run-up.

Comparisons are made again against the WCSPH method of Schwaiger & Higman [10] and the experimental results [2]. Both SPH methods used a similar particle distance to discretize the computational domain with d*x*=6 m and d*x*=7.5 m for the ISPH and WCSPH methods, respectively (referring to the full-scale problem). The corresponding particle count for the two methods is 10091 for the current work and 3677 for the WCSPH results of Schwaiger & Higman [10]. The difference in the particle count is partly due to the different particle spacing and initial volume of the landslide phase and partly because of the different boundary conditions used by the two methods (dummy particles for ISPH and particles with repulsive force for the WCSPH of [10]). This has, along with the increased computational cost of the divergence-free ISPH method and the addition of RANS turbulence modelling and non-Newtonian rheology (which also limits the time-step size), a significant effect on the computational cost of the proposed methodology; the ISPH simulation of the current work takes approximately 9 h for 60 s of physical time, compared with the 1 h of computational time reported by Schwaiger & Higman [10] on 3.4 and 3 GHz CPUs, respectively.

Figure 8 shows the comparisons for the wave run-up, the wave height at *x*=885 m and the landslide-entry profile between the experimental results of Fritz *et al.* [2], the WCSPH of Schwaiger & Higman [10] and the final ISPH results of the current work. It is clearly shown that with the proposed model a good overall agreement with the experiment of Fritz *et al.* [2] is achieved, having a close agreement for the wave run-up while retaining a good representation of the sediment-entry profile. Moreover, a significant overall improvement is noticed when compared with the results of Schwaiger & Higman [10], including both the good representation of the experimental landslide entry as well as an improved agreement in the wave run-up. Figure 8*b*, which corresponds to the wave height at position *x*=885 m, was the only measure found to under-perform against the results of Schwaiger & Higman [10], by failing to capture adequately the evolution of the wave in time period 14 s<*t*<21 s. Nevertheless, the agreement with the experiment is much closer for *t*>23 s which corresponds to the later stages of the wave propagation and the returning wave. A more thorough investigation of the physics at the impact, which may involve how the landslide particles are packed, should be considered to improve this apparent limitation.

Finally, figures 9 and 10 show the comparison of the final ISPH results and the experimental [2] for the wave run-up and landslide entry, respectively. Evidently, a close agreement with experimental findings is observed for the same time intervals. Noticeably, even the breaking of the wave before the run-up is captured in the ISPH results (figure 9*a*). Some minor differences are observed in the landslide shape during entry into the water phase, with some water particles being trapped below the landslide phase at the lower corner. Moreover, the water phase is found to be displaced by the landslide phase at a steeper angle compared with the experimental images, which may also partly explain the discrepancies observed in the wave evolution at *x*=885 m (figure 8*b*).

#### (iv) Effects of non-Newtonian, turbulence and saturation modelling

As shown in the previous paragraphs, the introduction of non-Newtonian and saturation models for the landslide phase as well as turbulence modelling for the water phase, along with true-to-the-experiment initial conditions, resulted in a close approximation of the experimental results for both the entry profile and the wave run-up. In this section, the effect of the introduction of the aforementioned models is briefly discussed. For further details, the reader is referred to Xenakis [31], where the effects of these models are discussed in detail.

#### (v) Introduction of non-Newtonian rheology

The introduction of non-Newtonian modelling is important in order to accurately capture the entry profile presented in the experiment of Fritz *et al.* [2], who used sediment to represent the landslide. When using the reconstructed pneumatic accelerator with Newtonian rheology for the landslide, it was found that the landslide entry-profile was over-predicted in elevation by an order of 10%, while having a misshapen profile. These findings, justify the shear-thinning behaviour expected by the sediment used in the experiment [2].

#### (vi) Turbulence modelling

The introduction of the *k*−*ϵ* turbulence model was considered appropriate since a high Reynolds number was measured with an order of 10^{9}. Moreover, with laminar regime in the water phase the run-up height was over-predicted by a factor of 3 when the pneumatic accelerator mechanism [52] was implemented. Indeed with the introduction of the *k*−*ϵ* turbulence model and the effect of eddy viscosity the run-up height was reduced by 20%. Nevertheless, this new run-up height was still more than twice higher compared with the experimental findings.

#### (vii) Saturation modelling

As explained in §2e and according to the findings of Mitarai & Nakanishi [44], a simple saturation model was introduced. The effect of this model is as described from Mitarai & Nakanishi [44] to transition the landslide phase towards a more shear-thickening behaviour, as it becomes wet. This transition and the increase of effective viscosity reduces the time for the landslide to come to a rest and therefore the rate with which the landslide is entering the water. Previous work by Walder [54] for non-breaking waves caused by subaerial landslides suggests that the wave amplitude can be predicted by the landslide flux per unit width. The resulting run-up height can be predicted by the wave-length and amplitude for a given depth (e.g. [55,56]). Although the aforementioned research cannot be used to predict the run-up height in the current case, since a breaking wave is involved, it can be seen that the reduction of the landslide flux induced by the saturation model can effectively reduce the height of the wave and therefore alter the run-up height. In summary, the saturation model presented herein both satisfies the rheological behaviour of the landslide from dry to wet state (as discussed also in [44]) and reduces the increased landslide flux, represented by the pneumatic accelerator mechanism in the experiments.

## 4. Conclusion

In the current work, the complex multi-phase environmental flow of the Lituya Bay landslide and tsunami has been examined. It was shown that previous computational methodologies failed to wholly capture the phenomena as presented in the benchmark experimental work of Fritz *et al.* [2]. Here, careful consideration of the initial condition as well as taking into account turbulence in the water phase and saturation of the landslide upon entry allows the accurate prediction of the benchmark findings. The method presented herein can be readily used to inform complex geotechnical applications which may involve non-Newtonian rheology, turbulence and saturation. Analysis of the results may then inform strategies that help prevent or prepare against the potentially devastating consequences that these geological events have on industrial activity and societal welfare.

## Data accessibility

Datasets and images are available from: http://dx.doi.org/10.6084/m9.figshare.4652572.

## Authors' contributions

A.M.X. built the computational model, carried out the data analysis, participated in the design of the study and drafted the manuscript; S.J.L. participated in the design of the study, helped draft the manuscript and participated in the coordination of the study; P.K.S. participated in the design of the study and helped draft the manuscript; B.D.R. participated in the design of the study, helped draft the manuscript and participated in the coordination of the study. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

A.M.X. has received financial support from the School of Mechanical Aerospace and Civil Engineering, University of Manchester.

## Acknowledgements

The authors would like to acknowledge the constructive comments of the members of the SPH group of the University of Manchester.

- Received September 7, 2016.
- Accepted February 20, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.