Open Access

Abstract. Tsunamis induced by rock slides plunging into fjords constitute a severe threat to local coastal communities. The rock slide impact may give rise to highly non-linear waves in the near field, and because the wave lengths are relatively short, frequency dispersion comes into play. Fjord systems are rugged with steep slopes, and modeling non-linear dispersive waves in this environment with simultaneous run-up is demanding. We have run an operational Boussinesq-type TVD (total variation diminishing) model using different run-up formulations. Two different tests are considered, inundation on steep slopes and propagation in a trapezoidal channel. In addition, a set of Lagrangian models serves as reference models. Demanding test cases with solitary waves with amplitudes ranging from 0.1 to 0.5 were applied, and slopes were ranging from 10 to 50°. Different run-up formulations yielded clearly different accuracy and stability, and only some provided similar accuracy as the reference models. The test cases revealed that the model was prone to instabilities for large non-linearity and fine resolution. Some of the instabilities were linked with false breaking during the first positive inundation, which was not observed for the reference models. None of the models were able to handle the bore forming during drawdown, however. The instabilities are linked to short-crested undulations on the grid scale, and appear on fine resolution during inundation. As a consequence, convergence was not always obtained. It is reason to believe that the instability may be a general problem for Boussinesq models in fjords.


Introduction
Rock slides and subaerial landslides are known triggers of large impulse generated tsunamis that may inundate coastal fjord communities.Although such tsunamis are rare, they may result in huge run-up in the vicinity of the landslide impact in the excess of those caused by earthquake tsunamis.Examples of rock slide induced tsunamis include the 1961 Lituya Bay event (Miller, 1960), the Lago Yanahuin (Plafker and Eyzagiurre, 1979), the 1783 Scilla landslide (Tinti and Guidoboni, 1988), and in 2007 a series of rock slides in the Aisén fjord in southern Chile caused tsunamis that were documented on video (Sepúlveda and Serey, 2009).In Norway, three major tsunamis struck the communities in Loen (1904Loen ( , 1936) ) and Tafjord (1934), causing altogether 175 fatalities (Jørstad, 1968;Harbitz et al., 1993).The tsunami hazard due to rock slides is significant in many communities in the western part of Norway (Blikra et al., 2005).A site that is considered particularly hazardous, is the unstable rock slope at Åknes in Storfjorden, where rock slide volumes of several million m 3 may impact the fjord.The Åknes rock slope is extensively monitored and displays relative movements of up to 20 cm yr −1 (Oppikofer et al., 2009).
Tsunamis induced by rock slides may involve a high degree of non-linearity in the generation process, including breaking, cavitation and strong turbulence, as clearly seen from scale experiments (Fritz et al., 2003;Saelevik et al., 2009;Mohammed and Fritz, 2012).In the far-field however, the dissipative terms are less important, and the wave propagation may be described by long-wave theory.Typically, both frequency dispersion and non-linearity are found to be Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
important.A possible modeling strategy could therefore be using a Navier-Stokes model in the generation area coupled to a depth averaged Boussinesq-type model for the wave propagation in the fjord system.The fjords are narrow, and are dominated by rugged steep slopes.Tsunamis inundate the coastlines as they propagate.At the same time the tsunami may exhibit breaking.Both of the latter effects should be properly accounted for in the propagation model, which obviously constitutes a challenge.
Due to the above mentioned capabilities, models such as FUNWAVE (Shi et al., 2012), and COULWAVE (Kim et al., 2009;Kim and Lynett, 2011) provide a state of the art modeling framework and have been popular for modeling dispersive tsunamis by landslides and volcanic flank collapses, (e.g., Lynett et al., 2003;Grilli and Watts, 2005;Geist et al., 2009;Abadie et al., 2012).On the other hand, such models were originally not developed for handling the violent flows due to tsunamis induced by rock slides or volcanic flank collapses involving strong non-linearity and simultaneous run-up along steep coastlines.Løvholt and Pedersen (2009) found that several Boussinesq formulations are prone to instability even in their linear formulation when subjected to steep bathymetric slopes.Stability issues are likely to be even more pronounced for highly non-linear waves in fjord systems where steep reliefs are present.Hence, there is a need to systematically test operational Boussinesq models to address their capabilities for simulating tsunamis under demanding conditions.
A fundamental requirement for any numerical model is the accuracy and convergence of the method.Convergence may be analyzed by means of grid refinement tests, demanding that the differences between the computed field variables vanishes as the grid lengths approach zero.A necessary and sufficient condition for convergence is stability (see e.g., LeVeque, 1992;Langtangen, 2003, for discussions).The first objective of this paper is to analyze convergence of operational models for a set of simple bench-mark problems for run-up on steep slopes.For this purpose, we will test convergence properties for one operational Boussinesq model (COULWAVE, Cornell University Long and Intermediate Wave Modeling Package), which exhibits the above mentioned properties.The second objective is to scrutinize how different run-up formulations deal with the steep slopes, combined with the strong non-linearity and dispersion.To this end, different run-up formulations are tested within the Boussinesq modeling framework.In addition to COULWAVE, we utilize a set of Lagrangian models that serve as reference.The Lagrangian models include both a boundary integral (BIM) full potential model as well as Boussinesq-type models.These models are mostly applicable to one direction of propagation, while some may be used in simple geometries with two horizontal dimensions, only.Hence, they are not an option for operational use, but are well suited for comparison because their good convergence properties are previously verified.The benchmark tests include inundation of a solitary wave on a slope.
The first part of this paper provides a review of run-up modeling using a depth-averaged framework, followed by a brief description the employed models.Section 3 presents a study of run-up on steep plane beaches using fully nonlinear non-hydrostatic models.Emphasis is put on model convergence and stability properties.This is followed by a study of wave propagation in a trapezoidal channel, which is presented in Sect. 4. Finally, the modeling capabilities of an operational Boussinesq-type model for landslide induced tsunamis in fjords is discussed.The breaking model employed in the operational model is presented in Appendix A. The Lagrangian reference models are described in references given below.In addition, particulars relevant for the present study are found in Appendix B.

Run-up modeling with depth integrated equations
The first theoretical treatment of non-linear run-up was presented by Carrier and Greenspan (1958) who transformed the non-linear shallow water equations on an inclined plane into the corresponding linear equations.They used the transformation to obtain a standing wave solution, as well as a run-up following from an initial condition.This technique has remained popular ever since and had a resurge after the formula for solitary wave run-up was published by Synolakis (1987).However, the technique is limited to hydrostatic equations and, save for a few exceptions (Kânoǧlu and Synolakis, 1998;Choi et al., 2008;Didenkulova and Pelinovsky, 2009) to plane slopes as well.Moreover, non-linear specifications of initial conditions are cumbersome due to the transformation technique, while the transformation back to the physical plane generally requires numerical integration.As a consequence only a few truly analytic solutions have been derived, and even the celebrated formula of Synolakis is an asymptotic approximation requiring a gentle slope compared to the incident wave length and becomes otherwise inaccurate (Pedersen, 2008a;Pedersen et al., 2013).In the present context, namely run-up of moderately short waves, we must instead rely on carefully obtained numerical solutions for comparison.
In Eulerian inundation models computational cells must be redefined as wet or dry according to the shoreline motion.This often requires a special treatment of near-shore points and may involve extrapolation of field quantities to newly flooded cells or fictitious grid points.Sielecki and Wurtele (1970) published the, maybe, first proper attempt on such modeling.Later, Hibberd and Peregrine (1979) employed the Lax-Wendroff method, combined with a multi-step scheme for advancing the shoreline for bores as well as non-breaking waves.Today, most standard models for tsunami applications or coastal engineering come with some form of inundation facility.The TUNAMI (Imamura, 1996) model is based on the NLSW equation and employs a stair-step procedure in the sense that the depth is regarded uniform in each cell and that a dry cell becomes flooded when its shelf is overtopped by the fluid elevation in a neighboring cell.MOST (Titov and Synolakis, 1998) is another widespread NLSW solver invoking a split step method with alternating directional application of characteristics.The shoreline is traced by an auxiliary grid point, while values in newly flooded cells are projected from the neighboring wet nodes.In NLSW models based on Riemann solvers, yielding TVD schemes, the shoreline may be implemented as a special Riemann solution (propagation into vacuum) combined with requirement on minimum flow depth (see, for instance, LeVeque and George, 2008).Bellotti and Brocchini (2001) invoked the TVD scheme and the Riemann shoreline technique in a Boussinesq framework.This has recently been adopted in the standard Boussinesq models, such as FUNWAVE and COULWAVE.The older versions of the FUNWAVE (Kennedy et al., 2000) model for Boussinesq-type equations come with a particular wet-slot treatment of the beach, which then displays properties akin to a porous medium.In recent descendants of the FUNWAVE models (Shi et al., 2012) the authors have resorted to more standard techniques (Tonelli and Petti, 2012).Another standard Boussinesq model (Lynett et al., 2002) employed extensive onshore extrapolation of field variables to reduce the need for special treatment of shoreline points.Generally, in Boussinesq-type models with moving shorelines the dispersion term is partly or fully deleted in the vicinity of the shoreline.Herein, we employ a descendant of the COULWAVE model (Kim et al., 2009;Kim and Lynett, 2011) referred to as the operational model below, combined with different techniques for run-up.Presumably, the obtained results will be relevant also for the application of these techniques with other basic numerical models.
A generally simple and robust way to deal with a moving shoreline is to transform the basic equations to a coordinate system that deforms accordingly.This may be obtained by applying Lagrangian coordinates (Goto, 1979;Goto and Shuto, 1983;Pedersen and Gjevik, 1983;Zelt and Raichlen, 1990;Johnsgard, 1999) or a more flexible ALE (arbitrary Lagrangian-Eulerian) description ( Özkan-Haller and Kirby, 1997;Prasad and Svendsen, 2003).For such models rapid convergence may often be obtained, but they are limited to moderately complex geometries and are thus not suitable for tsunami run-up, for instance.On the other hand, since they are accurate and offer some freedom concerning geometry and physical description they may serve excellently for benchmarking more general operational models.Unfortunately this is not much exploited.Due to tradition, maybe, authors tend to validate their models by experiments, which are hampered by scaling effects and often issues concerning the definition of the problem, or by the so-called analytical, but certainly limited, solutions obtained from the hodograph technique or the few other analytic solutions which do exist, such as oscillations in parabolic basins (Thacker, 1981) and dam-break (Stoker, 1957).The widely used models of today are generally tested on a set of such problems with good results.However, as the range of the available analytical solutions is quite limited, this offers no guarantee for good performance for even moderately more complex problems.This will be demonstrated in the present manuscript where we focus on a run-up problem, involving strong nonlinearity, dispersion and steep slopes, and propagation in a channel of non-rectangular cross section.Herein, we do use Lagrangian techniques for comparison.One set of models is based on Boussinesq-type equations, with some diversity with respect to non-linearity and dispersion properties (Jensen et al., 2003;Pedersen, 2008bPedersen, , 2011;;Pedersen et al., 2013).Since these models possess very accurate shoreline tracing, they easily yield numerical solutions very close to convergence.The diversity between the different Boussinesq models yields a (generally narrow) range of run-up heights, which indicates what may be expected for Boussinesq solutions in general.In addition we employ a boundary integral technique for full potential theory (Pedersen, 2008b;Pedersen et al., 2013), without any approximations with respect to wavelength or non-linearity.While the Lagrangian long-wave models are mainly used for assessing the numerical performance of operational models, the boundary integral method provides a check on the physical validity as well.

Employed Boussinesq models
We introduce a Cartesian coordinate system with horizontal axes, ox and oy in the undisturbed water level and an oz axis pointing vertically upward.The equilibrium depth is denoted by h, the surface elevation by η and the velocity components by u and v in the x and y directions, respectively.We identify a typical depth, d, a typical wavelength, L, and an amplitude factor, , which corresponds to a characteristic value of η/d.Different long-wave equations can be obtained through perturbation expansions in µ ≡ d/L and .They may then www.nonlin-processes-geophys.net/20/379/2013/ Nonlin.Processes Geophys., 20, 379-395, 2013 be classified according to which orders these parameters are retained in the equations, when the equations are scaled such that the leading order is unity.The residual (error) terms of the standard Boussinesq equations, such as solved in the early Boussinesq models (Peregrine, 1967), are O( µ 2 , µ 4 ).The primary unknowns then were the surface elevation and the vertically averaged horizontal velocity.Several other formulations with different choices of primary unknowns do exist, of which that of Nwogu (1993) has become widely used.
In this formulation the velocity at a chosen depth z α is used as a primary unknown.With the optimal choice z α = −0.531himproved linear dispersion properties are obtained (good for, say, wavelengths down to 2h).Wei et al. (1995)  ).These equations have later been corrected and generalized to include multiple layers and turbulent shear effects (Lynett and Liu, 2004;Kim et al., 2009;Kim and Lynett, 2011).
Herein we employ different varieties of an operational model (COULWAVE) and a set of reference models that have Lagrangian shoreline tracking as common feature.The particulars of the models will be explained subsequently, but it is convenient to introduce a brief definition with a numbering already at this stage.The numbering will be used for reference later in the text and in figures: e. boundary integral model based on full potential theory.

Eulerian operational model
The COULWAVE model was first developed as a means to investigate waves generated by submarine landslides, and was numerically very similar to the initial versions of the FUNWAVE model.Recently, the numerical scheme has been changed to utilize a finite-volume (FV) method for the Boussinesq equations in conservative (flux) form (Kim et al., 2009).Various turbulence and rotational effects have also been included (e.g., Kim and Lynett, 2011), but these features are not utilized here.With the change to the FV scheme, the moving boundary approach was also modified, discarding the earlier "extrapolation" technique (Lynett et al., 2002).While the extrapolation method proved generally stable and accurate, it was unable to handle complex flow convergence and flow re-entry.Following existing FV moving shoreline approaches for flux-form equations, the wet/dry boundary can be accommodated numerically through special treatment of the fluxes in the boundary cells.Consider a situation where cell i is wet and cell i + 1 is dry; here, the fluxes at the interface, i + 1/2, require special calculation.The flux terms for a one-horizontal-dimension configuration are the mass flux (H U ) i+1/2 and the momentum flux (H U 2 ) i+1/2 , where H = η +h is the total water depth and U is the depth-averaged velocity.There will be three different approaches for estimating these terms presented in this paper, given below.Note that for all three run-up approaches, the boundary cell flux terms H and U are handled independently, such that, for example, the mass flux becomes H i+1/2 U i+1/2 .
(a) Centered method.In this case, H i+1/2 = 0.5(H i + H i+1 ) and U i+1/2 = U i .Here, for all slope and flow depth configurations, η information from the first dry cell is included in the boundary flux calculation, and conceptually the flow depth and depth profile vary linearly across the wet and dry cells.
(b) Step method.For this approach, the conventional stairstepped schematic of the bottom is used.For interface fluxes to be non-zero, η i must be greater than −h i+1 .When this condition is satisfied, H i+1/2 = η i + h i+1 and U i+1/2 = U i .This is a low-order approach, and is the method described in Lynett et al. (2010, with details on non-simple wet/dry cell configurations and other details of the moving boundary approach).The depth profile is imagined to step up (or down) vertically at x i+1/2 from h i to h i+1 , and the free surface gradient in this region is zero.
(c) Hybrid method.The concept here is very similar to the stepped flux approach in (b), except that the depth profile between the center of cells i and i + 1 is considered to vary linearly, as in (a).For interface fluxes to be non-zero, η i must be greater than For all three approaches, cell velocities at non-boundary dry cells are set to zero, and the free surface gradient at the last wet cell is evaluated with a low-order directional difference away from the dry cell, i.e., (∂η/∂x) i = (η i −η i−1 )/dx.All dispersive terms at the last wet cell are neglected.

Lagrangian run-up models
The Lagrangian models are particularly designed to deal with run-up on sloping beaches for simple geometries and enable a shoreline description of high accuracy.A fully non-linear extension of the standard formulation (Pedersen and Gjevik, 1983;Jensen et al., 2003;Pedersen, 2008b)  ), which makes the equations comparable to those of Wei et al. (1995).Different versions, valid to this order, are obtained through a free parameter κ.Linear dispersion properties are improved by adding an O(µ 4 ) term to the momentum equation and then optimizing the linear dispersion properties in the same manner as Nwogu (1993).More details are given in Appendix B, including a precise definition of the model variants 2a, 2b, and 2d (see above).
For generality we also use a combined Boussinesq/NLSW model, named 2c, as described by Pedersen (2011).In finite depth this model employs standard Boussinesq equations discretized on an Eulerian grid, while a small Lagrangian grid is used near shore.The different Lagrangian models are used to obtain converged numerical run-ups and, by means of their diversity, to indicate a range for Boussinesq results.As a check on the accuracy of the long-wave approximation as such we also compare with the results of a boundary integral model for full potential theory (Pedersen, 2008b), named 2e.

Plane wave run-up on steep slopes
The model's capabilities of reproducing run-up on steep slopes are investigated by examining solitary waves on inclined planes.The incident waves are specified by the procedures described below, while the bathymetry consists of a flat bottom joined smoothly to the inclined plane by means of a spline function.More precisely, if the intersection between the flat bottom and the inclined plane is located at x = x k the depth in the interval − < x − x k < is given by a polynomial of fifth degree which yields continuity for h, dh/dx and d 2 h/dx 2 .The spline is invoked to avoid effects of a vertex in the bottom profile.Even though such effects may be important in their own right, our intention with the present test is to study model performance for run-up on steep slopes in a context as simple as possible.
We choose d as the equilibrium depth on the flat bottom, giving a dimensionless depth smaller than or equal to unity.The employed slope inclinations are θ = 10 • , 15 • , 20 • , . . ., 50 • and is d•cot θ/10.A solitary wave is characterized by A/d, its amplitude to depth ratio.Herein we employ solitary waves with amplitudes ranging from A/d = 0.05 to A/d = 0.5.Synolakis (1987) combined a linear treatment of an incident solitary wave with a non-linear shallow water theory on the slope to obtain a celebrated formula for the relative runup height, R/A, as function of A/d and θ, as well as a breaking criterion.However, this approach is not accurate for steep slopes such as studied herein (see, for example, Pedersen, 2008a).Hence, we must rely on carefully performed simulations and comparison between models.Grilli et al. (1997) used a boundary integral method for solitary waves incident on a plane.No maximum run-up heights were reported, but a criterion for breaking during run-up was fitted to the com-putational data: where C = 25.7 indicates clear overturning and C = 16.9 the lower limit of an intermediate regime where very steep fronts were formed.For θ = 10 • Eq. (1) yields A/d ≥ 0.52 for C = 16.9, while it predicts that no incident solitary wave breaks during run-up for angles larger than 12.5 • .Hence, according to this criterion none of our solitary waves should break during run-up, even though they generally break during drawdown.

Model setup and incident waves
Different Boussinesq equations inherit different solitary wave solutions.Only some of the formulations described herein possess known closed form solitary wave solutions, namely the non-linear Lagrangian forms, without optimized dispersion (see Appendix B).For full potential theory we have the numerical solution of Tanaka (1986), while a perturbation solution, namely the fourth-order solution of Fenton (1972), may be applied for small amplitudes (A/d < 0.1).A numerical solution is also employed for the standard Boussinesq equation, employed in model 2c (see Pedersen, 1988).
Then, for Boussinesq models with improved dispersion properties we have only access to approximate analytic solitary wave solutions.For the fully non-linear Lagrangian model with improved dispersion properties this solution is given in Appendix B, while the approximate solution for the equations of Wei et al. (1995) is found in Wei and Kirby (1995).
For larger amplitudes the deviation between the approximate and the exact, undescribed, solution becomes significant.
When an approximate solitary wave solution is inserted in a model its height and shape will gradually adapt, while a tiny wave train is shed at the rear of the leading pulse.If the model is without damping the amplitude and shape may approach a stationary state corresponding to an exact solution consistent with the model (see Pedersen, 1991).If there is a weak damping the model will produce a slowly attenuated wave with properties that otherwise are very close to those of a solitary wave solution of the underlying equations.
For each model we employ initial conditions which represent the exact solitary wave solution for the amplitude given, if such a solution is available.For the models which lack an exact solution we insert the approximate solution, with amplitude A 0 , as initial conditions in a very fine grid, and propagate the solution a distance L c over constant depth.We then obtain a wave with amplitude A L and a shape which is closer to a perfect solitary wave owing to the equation set in use.When this procedure is repeated for a number of different A 0 we obtain a table of A L as function of A 0 , which we may invert through interpolation to obtain A 0 (A L ).In the run-up simulations we then employ initial conditions corresponding to an approximate solution with amplitude A 0 at a distance L c in front of the start of the beach.When the wave reaches the toe of the beach it is then close to a solitary wave of height A. We have chosen L c = 30d.Naturally, the procedure is not exact and even though the amplitudes are the same at the start of the shoaling process for all models, the wave shapes depend on the particulars of the individual models.
The relation between A 0 /d and A L /d for a range of amplitudes is depicted in Fig. 1.For the COULWAVE formulation the solution of Wei et al. (1995) for the initial conditions is employed.The grid resolutions range from x/d = 0.06 to x/d = 0.02, while the Courant number is 0.1 and 1 in the COULWAVE model and the Lagrangian model, respectively.In the run-up simulations below, we use the inverse relation to provide a look-up table for the initiation of the wave field.Examples of different slopes and initial waves are found in Fig. 2. numbers, respectively.The run-up R is normalized by the initial surface elevation at the slope start (R/A).For the Lagrangian models, results are displayed only for the finest resolution.It is noted that convergence tests for the various Lagrangian models generally showed high accuracy; the reported maximum run-up compared with those from simulations using 1-2 times lower grid resolution generally provided errors of less than 0.1 %, usually even 1-2 orders of magnitude smaller.Exceptions are the combinations of high amplitude with 10 • slope and some simulations at the steepest inclinations.From Fig. 3, we also observe that there is apparently a very little spread in the run-up between the different Lagrangian reference models.For the operational Boussinesq model, there is significantly more scatter, and for certain parameter combinations their results differ significantly from those of the reference models.The largest deviations are due to model instability before the maximum run-up was reached, observed for the largest amplitudes and for the smallest slope angles in particular.A coarser spatial grid resolution (Fig. 3), and smaller Courant number (Fig. 4) provide more stable results.The accuracy of R/A using model 1a relative to reference model 2d is displayed as a function of the different slopes and amplitudes in Fig. 5. Three different grid resolutions and a Courant number of 0.1 are employed.The relative spread in the R/A ratios for the reference models are shown for comparison.Figure 5 also displays for which parameter combinations instability is reached.Generally, the coarse resolution provides stable results but yields much lower accuracy than the reference models.The operational model is most accurate for small values of A and large slope.Furthermore, the accuracy is improved with increasing grid resolution but here the model also becomes more unstable.It is also more prone to instability for increasing A and decreasing θ.It is noted that for the highest waves the finest resolution was almost 40 grid points per depth.This is considered extremely fine, and most likely not feasible to implement in a real application in 2HD (two horizontal dimensions).The results displayed in Figs.3-5 are obtained from simulations using a flux limiter.However, limiters are mainly employed to amend stability problems and may affect the accuracy.Figure 6, compares R/A ratios both with and without use of limiters.The simulations without limiter provide a few percent higher run-up on some occasions, but are more often indistinguishable from the cases where the limiters are used.Figure 6 alone would indicate that the simulations are not very sensitive to the use of limiters.For the run-up of high waves on steep slopes examined here, the flux limiter did not improve the stability significantly.However, the use of limiters provides more stable results in other situations, for instance during drawdown.

Convergence tests using different run-up formulations
Further investigations using the operational Boussinesq models 1b and 1c were undertaken to examine the effect of the numerical run-up formulation.Simulations were conducted with different slope angles and values of A/d.R/A values for the finest grid resolutions are depicted in Fig. 7, results using method 1a and the Lagrangian models are retained.More stable results and closer agreement with the reference model are obtained.The deviation from the reference solution increases with larger initial amplitude A/d.Method 1c provided a particularly good match.However, from the authors' experience method 1c is less robust than 1a and 1b in 2HD with more complex geometries.

Wave evolution and spurious effects on the slope
From the convergence tests it was found that instability occurred during the first positive inundation.The simulated surface elevations prior to instability using COULWAVE using method 1a is depicted in Fig. 8 for A/d = 0.3 with different slopes and grid resolutions.We have observed that instability occurs for the gentlest slopes tested combined with high grid resolutions.Instability before maximum run-up is displayed in the right panels in Fig. 8, both for high spatial resolution x/d = 0.03 and gentle slopes (θ = 10 − 15 • ).
Running the simulation at a coarser resolution of x/d = 0.06 or at a steeper slope (θ = 20 • ) yielded instability during withdrawal.We observe that the breaking facility is invoked already during run-up, although the waves should be non-breaking (see discussion below and Eq. 1).Moreover, the breaking facility is not effective during withdrawal in one occasion, due to the shore-facing bore being stationary.tested.Trying other breaking facilities (Kennedy et al., 2000) caused instability to occur at coarser resolutions and at earlier stages of propagation.reflections and scattering may be more pronounced.The fjords are narrow, often only a few kilometers wide, and the wave propagation is expected to be heavily affected by the simultaneous run-up along the steep coastline.As these steep parts of the fjords are often not inhabited, a precise simulation of run-up itself is not our primary interest.Rather, we are interested in the accuracy of the model for describing the distant wave propagation.
The trapezoidal channel constitutes a simplified fjord topography, or half the topography of a symmetric fjord.Lynett et al. (2002) studied solitary wave propagation in a trapezoidal channel with steep slopes using COULWAVE.Simulated surface elevations and inundation compared favorably against laboratory experiments conducted by Peregrine (1969), although with a minor underprediction of the run-up.
In the experiments of Peregrine (1969) slopes were 1 : 1 and the width of the channel was restricted to 1.5d.Moreover, solitary wave amplitudes were restricted to A/d < 0.2.In the following, we investigate the solitary wave propagation in a trapezoidal channel that resembles a typical fjord.The channel has a maximum water depth of 150 m, the constant depth part a width of 2 km, and the side slopes are 13 • .The channel is depicted in Fig. 9. Figure 9 also shows the propagation of a solitary wave with A = 27 m at different time steps using model 1b.Due to symmetry, only one half of the channel is included in the simulation.Furthermore, the figure shows that even for this relatively simple configuration, the initial solitary wave shape is rapidly distorted.Along the sides, the waves are retarded and the wave is refracted, producing a second wave.Hence, the run-up process affect the propagation along the channel.This process repeats itself causing rather complicated wave dynamics, which are elaborated in Lynett et al. (2002).
A grid refinement test was conducted for the geometry in Fig. 9.In these simulations, no bottom friction was assumed, but the breaking facility was invoked.Three different solitary wave amplitudes A/d = 0.18, A/d = 0.3, and A/d = 0.5 were considered.Transects of the surface elevation were compared at the center of the channel, indicated by the white dashed line in Fig. 9. Stable results for model 1b for non-linearities A/d = 0.18 and A/d = 0.3 are depicted in Fig. 10 at different time steps and grid resolutions.Similar simulations conducted for model 1a provided stable results for the same choice of model parameters, but the results were much less accurate as also indicated by the run-up simulations above.For the finest grid resolutions, the deviation of the maximum surface elevation is less than 5 %.For this degree of non-linearity, the model seems to converge properly and be well suited.For the largest non-linearity however, the simulations broke down except for relatively coarse grid resolutions.Courant numbers of 0.25 and 0.1 were used, and the smallest time steps provided improved stability.This was particularly important for the highest non-linearities.It is expected that instabilities linked to the inundation are the cause of the model break down.However, unstable wave shapes similar to those shown in Fig. 8 were not found.It is stressed that an amplitude of A/d = 0.3 is considered truly non-linear, and that the model seems to be able to deal properly with this in the present geometry.

Concluding remarks
Instabilities are revealed for highly non-linear solitary waves in steep slope geometries for an operational Boussinesq model.A consequence of the instability is that accurate modeling of run-up may not be feasible when the wave amplitude is too large.Instabilities arise most easily during drawdown, and for run-up under more demanding conditions.They are triggered at fine spatial grid resolutions, and therefore the simulations fail to converge.Inclusion of breaking terms or limiters, or reduction of the time steps, only delayed the onset of the instability.The stationary bore formed during drawdown becomes particularly challenging, as the breaking criterion is not effective.In other situations, false breaking occurs prior to instability.Replacing the centered run-up method with either a stepped or a hybrid method provided much increased accuracy.It is noted that the hybrid approach is not a feasible strategy as this method is more prone to instability in real geometries.Simulated wave propagation and simultaneous inundation along the steep side slopes of the trapezoidal channel show that the propagation is affected by the inundation, and that wave energy from the leading wave train is propagated backwards and dissipated due to run-up.When the non-linearity and grid resolutions exceed a certain level, instabilities are triggered.For A/d = 0.3 the model gives stable results, but instability is observed for A/d = 0.5.For a tsunami induced by a large subaerial landslide in a fjord, A/d will normally exceed 0.3 in the near field.However, wave amplitudes will attenuate relatively fast to A/d < 0.3.A/d > 0.3 therefore represents a relatively demanding test for the far-field propa- non-linearity occurs, to amend instabilities.Coupling of a primitive near-field model to a Boussinesq model for the farfield propagation may be the best option.One example of such a coupling was demonstrated by Løvholt et al. (2008) for the modeling of a potential tsunami from the La Palma Island.As this model was one-way only, we suggest that future models for tsunamis induced by subaerial landslides should be two-way coupled.
The instabilities occur also in tests with similar models other than COULWAVE (results not shown).This suggest that the instability is a general problem that may occur in any operational Boussinesq model attempting to include run-up.Still, the different run-up formulations investigated provide different stability properties and accuracy, but none of them removed the instability altogether.It is furthermore stressed that the reference models, which are based on a different mathematical formulation, provide more coherent results.However, operational models cannot be based on these formulations.So far the operational models have only been tested in idealized geometries, but experience tells us that instabilities are triggered more easily when subject to complex geometries.Here, instabilities may arise differently compared to the idealized cases.Hence it is important to further simulate the tsunami propagation for a real fjord.

Appendix A The transport based breaking criterion
The primary motivation for developing another breaking scheme for the Boussinesq equations was an inadequacy for a wide range of simulations with the existing schemes.The primary detraction of these existing schemes is that the breaking dissipation, provided through an eddy viscosity, is entirely local; while the breaking event age is incorporated into various schemes, there is no permitted advection or diffusion of breaking induced turbulence.This essentially means that the viscous memory of the breaking dissipation is very weak.This does not seem reasonable in the surf and swash zone, and typically other dissipation mechanisms (e.g., subgrid or bottom friction) must be manipulated to provide good results in these regions.
Thus, here an attempt is made to model the breaking induced turbulence with a set of transport equations.The reader should keep in mind that detailed small-scale breakingdriven physics, which by definition are disregarded through the fundamental assumptions of the Boussinesq-type derivation, should not be properly provided by this model.This is, as are all other breaking models in depth-integrated equations, a practical and ad hoc submodel added to the governing equations for engineering purposes.
The transport-based breaking model aims to provide an eddy viscosity, and that eddy viscosity is added to the governing momentum equations through mixing terms.The for-mulation of the mixing terms is the same as that presented in Kennedy et al. (2000).The eddy viscosity is calculated through a local transport equation: where ν is the breaking-induced eddy viscosity, H is the total water depth, c = √ gH is the local non-linear long-wave speed, and ν source is the additional "source" of eddy viscosity arising from local breaking.Note that this equation is a local equation, meaning there is no apparent transport of any information.In fact the transport effects arise only from the ν source term.This term is given by where B provides a dimensionless measure of the local breaking intensity and η is the water surface elevation.The transport of information in this model is through the B term only.A simple advection equation with added source and diffusion term is used: where the p provides the B source, and will be described in a moment.In the above equation, the first two terms on the right hand side are the transport terms, which are upwind differenced using the local flow velocity to determine the upwinding direction.The formulation of the source term p relies largely on existing breaking threshold studies, manipulated slightly here to allow for smoother breaking initiation and cessation.While in existing models, breaking turns on and off based on a binary threshold, here we use tanh functions to control the "start" and "end" of a breaking event.For our breaking parameter, we will use the temporal change in the free surface elevation:  ) This equation can be interpreted as follows.In a region where breaking has not yet started (B = 0), the threshold for dissipation to growth will be Q = 0.65, while in regions where breaking is ongoing (B large), the threshold for breaking to end is Q = 0.40.In existing studies, this breakingevent variation of the threshold is typically dependent on the breaking age.While this is somewhat implicit in the formulation used here, more precisely our breaking threshold is dependent on the intensity of the breaking dissipation.We have found this to provide more reasonable wave-shape predictions across a range of different setups.Finally, the B-source term, p, is formulated as which again uses the tanh functions to provide a smooth forcing.We can see here that when Q R, then p goes to zero, and there is no ongoing local breaking event, while when Q R, then p goes to one, and the breaking event is ongoing and strong.The above set of equations comprises the complete formulation for this scheme.Note that the various integer coefficients in the equations are all tuned parameters.Included in this description are the comparisons for five regular wave cases, as presented in Kennedy et al. (2000) (original data from Hansen and Svendsen, 1979).Numerical results are shown in Figs.A.1-A.5.The comparisons are as good as or better than existing comparisons in the literature.Note the grid convergence tests for these trials, showing the model's ability to stably run breaking simulations with extremely high resolution.

Appendix B Lagrangian long-wave model
In this appendix we adopt the standard scaling for evolution of long-wave equations.Dimensionless coordinates (marked by ˆ) are defined according to where g is the acceleration of gravity and is as defined in Sect.2.2.Since this scaling will be used systematically while a fully non-linear momentum equation may be expressed: (1 +P 2 + γ P 3 + O(µ 4 ). (B3) The implicit parts P 1 and P 2 are additional non-hydrostatic terms When κ is set to unity and P 3 omitted, the set (B2) and (B3) equals the Boussinesq formulation of Jensen et al. (2003).
Collection of the κ terms on the left-hand side yields the net contribution κµ 2 h ∂η ∂x ( Du Dt + ∂η ∂x ), which is O(µ 4 ) since the two terms within the parenthesis cancel out to O(µ 2 ) (leading order balance of the momentum equation).factor κ may be chosen freely without violating the formal order of the equations, but it may change the properties qualitatively for steep waves.For κ = 1 we experience a singularity when the contact angle at the beach reaches 90 • , which is consistent with full potential theory as explained in Jensen et al. (2003), while no such singularity arises for κ = 0.However, this distinction is not relevant for the tests run herein and we employ κ = 0 and κ = 1 only to find a range of runup values for Boussinesq-type equations.
The term P 3 is This is a non-linear counterpart to the correction term suggested by Madsen and Sørensen (1992).For linear waves over constant depth the choice γ = 1 3 + z α h + z 2 α 2h 2 reproduces the dispersion relation of Nwogu (1993).The optimal choice z α = −0.531then corresponds to γ = −0.057.Madsen and Sørensen (1992) proposed instead the choice γ = − 1 15 , which reproduces the first three terms of the Taylor expansion of the dispersion relation.
The transformation to the Lagrangian coordinate, a, is then defined according to an averaged particle, moving with u: ∂x(a, t) ∂t = u, x(a, 0) = a . (B4) The continuity equation is then integrated to H ∂x ∂a = H (a, 0) ≡ H 0 , (B5) which expresses volume conservation in a material fluid column.The momentum equation may be re-casted into the form This equation is expressed partly in conservative form, which is convenient rather than necessary for non-breaking waves.
The bottom pressure term H 0 {h } must be given a representation that matches the preceding term on the right-hand side at equilibrium.Some further details on the discretization are found in Jensen et al. (2003).
The definition of models 2a-c in Sect.2.2 may now be expressed.
A solitary wave is a wave of permanent shape and constant wave celerity, c.
with Lagrangian shoreline tracking: a. weakly non-linear Boussinesq with standard dispersion; b. fully non-linear Boussinesq with standard dispersion; c. domain decomposition model (lagrangian run-up model); d.Boussinesq model, fully non-linear with optimized dispersion;

Fig. 1 .
Fig. 1.Initial solitary wave amplitude as a function of the amplitude after a propagation distance of x/d = 30.The very close agreement between two of the curves is presumably a coincidence.
Numerical simulations are made for a range of initial conditions at the base of the slope, ranging from A/d = 0.05 to A/d = 0.5, with steps of A/d = 0.05.We first present results for model 1a.Grid resolutions of n = 100, 200, and 400 points per wavelength, and Courant numbers of 0.1 and 0.5 were used for the spatial and temporal discretization, respectively.The transport-based model outlined in Appendix A is used to invoke dissipation due to wave breaking.Corresponding simulations using different Lagrangian models are conducted for comparison.The maximum run-up were extracted from the simulations, and are displayed as a function of the slope in Figs.3-4, indexed by their different grid resolutions and Courant Examples of the model setup, for three different initial wave configurations and bathymetric slopes, respectively.

Fig. 3 .
Fig. 3. Run-up ratio for solitary wave run-up using model 1a for different slopes, initial amplitudes, and grid resolutions.Here, a Courant number of Cr= 0.1 is employed.Results with Lagrangian models are included for comparison.Upper panel, A/d = 0.1, mid panel A/d = 0.3, lower panel A/d = 0.5.The drop in the run-up ratio for the largest amplitudes in combination with gentlest slope is due to model breakdown.

Fig. 4 .Fig. 5 .
Fig. 4. Run-up ratio for solitary wave run-up using model 1a at different slopes, initial amplitudes, and time steps.Here, a grid resolution n = 200 points per wavelength is employed.Upper panel, A/d = 0.1, mid panel A/d = 0.3, lower panel A/d = 0.5.The drop in the run-up ratio for the largest amplitudes in combination with gentlest slope is due to model breakdown.

Fig. 6 .
Fig. 6.Run-up ratio for solitary wave run-up using model 1a at different slopes.Results both with and without use of limiters are displayed.The Courant number is 0.1 and the spatial resolutions resemble those in figure 4 with n = 200 points per wavelength.Upper panel, A/d = 0.1, mid panel A/d = 0.3, lower panel A/d = 0.5.The drop in the run-up ratio for the largest amplitudes in combination with gentlest slope is due to model breakdown.

Fig. 7 .Fig. 8 .
Fig. 7. Run-up ratio for solitary wave run-up in COULWAVE for different slopes and initial amplitudes.Upper panel, A/d = 0.1, mid panel A/d = 0.3, lower panel A/d = 0.5.Three different run-up methods are employed.Here, a Courant number of Cr= 0.1 is employed.Dropout or underprediction of the R/A values for the gentlest slopes indicate instability.Results with Lagrangian models are included for comparison.

Fig. 9 .
Fig. 9. Upper panel, employed bathymetry for the trapezoidal channel.Due to symmetry, only one half of the channel is depicted.The color bar gives the total depth in m, negative numbers indicating points below still water level and positive values indicate land.The dashed white line indicated the location for the evaluation of the transects.Lower panel, surface elevation in m at different times for A = 27 m.Here, the color bar indicates the water elevation in m.
Fig. A1.Computed and measured wave heights and setup for the Hansen and Svendsen case 031041 (kh o = 0.369, H / h o = 0.12).
Fig. A2.Computed and measured wave heights and setup for the Hansen and Svendsen case 041041 (kh o = 0.501, H / h o = 0.11).
Fig. A3.Computed and measured wave heights and setup for the Hansen and Svendsen case 051041 (kh o = 0.641, H / h o = 0.11).