How does multiscale modelling and inclusion of realistic palaeobathymetry affect numerical simulation of the Storegga Slide tsunami?

The 8.15 ka Storegga submarine slide was a large ( 3000 km), tsunamigenic slide off the coast of Norway. The resulting tsunami had run-up heights of around 10–20 m on the Norwegian coast, over 12 m in Shetland, 3–6 m on the Scottish mainland coast and reached as far as Greenland. Accurate numerical simulations of Storegga require high spatial resolution near the coasts, particularly near tsunami run-up observations, and also in the slide region. However, as the computational domain must span the whole of the Norwegian-Greenland sea, employing uniformly high spatial resolution is computationally prohibitive. To overcome this problem, we present a multiscale numerical model of the Storegga slide-generated tsunami where spatial resolution varies from 500 m to 50 km across the entire NorwegianGreenland sea domain to optimally resolve the slide region, important coastlines and bathymetric changes. We compare results from our multiscale model to previous results using constant-resolution models and show that accounting for changes in bathymetry since 8.15 ka, neglected in previous numerical studies of the Storegga slide-tsunami, improves the agreement between the model and inferred runup heights in specific locations, especially in the Shetlands, where maximum run-up height increased from 8 m (modern bathymetry) to 13 m (palaeobathymetry). By tracking the Storegga tsunami as far south as the southern North sea, we also found that wave heights were high enough to inundate Doggerland, an island in the southern North Sea prior to sea level rise over the last 8 ka. 2014 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/3.0/).


Introduction
Around 8150 years ago the Storegga submarine slide generated a large tsunami that spread across the Norwegian-Greenland sea (Haflidason et al., 2005;Bondevik et al., 2005;Løvholt et al., 2005). The submarine slide had a volume of between 2400 and 3200 km 3 , affecting a region of 95,000 km 2 , making it one of the world's largest exposed submarine slides (Haflidason et al., 2005). The volume of material within the Storegga Slide is around 300 times the modern global annual sediment flux from rivers to the oceans. The Storegga slide is bigger than Scotland, and its headwall extends for $300 km. It dwarves even the largest slide yet found on land.
Many tsunami deposits from the Storegga slide-generated wave have been found across the region, including Scotland (Smith et al., 2004;Tooley and Smith, 2005;Dawson and Smith, 2000;Long et al., 1989;Dawson et al., 1988) northern England (Boomer et al., 2007), Norway (Svendsen and Mangerud, 1990;Bondevik, 2003;Vasskog et al., 2013) Faroe Islands (Grauert et al., 2001), and Greenland (Wagner et al., 2007). Run-up heights are estimated to be over 20 m in some locations, particularly where the tsunami wave propagated large distances along Norwegian fjords (Vasskog et al., 2013). The Storegga slide is the only large slide-tsunami that has been mapped out in such detail and over such a large area. This makes it an ideal case-study to examine basin-scale tsunamigenic slides. Numerical simulations of the slide and the subsequent tsunami (Harbitz, 1992;Bondevik et al., 2005) have shown how the wave propagated and are in reasonable agreement with run-up heights inferred from geological observations. However, previous models have been limited by two important technical constraints. First, they used relatively low spatial resolution along coastlines due to the large region simulated. This means that wave propagation along the complex Norwegian coast, for example, may not be properly simulated. Second, all previous studies used modern http bathymetry, as opposed to the inferred bathymetry from 8150 years ago, which has likely changed by tens of metres as a result of non-uniform isostatic relative sea-level changes.
Numerical simulations are a useful tool for studying tsunamis. A number of previous studies have used numerical models to study land-and submarine-slide generated tsunamis (e.g. Abadie et al., 2012;Assier-Rzadkieaicz et al., 2000). They allow some quantification of the hazard posed by such events, which is uncertain . A number of these studies have used nested models (multiple, coupled models with different spatial resolutions, using one or more codes) to simultaneously simulate both the large region and local details (Allgeyer et al., 2013;Kirby et al., 2013;Horsburgh et al., 2008). In particular, Bondevik et al. (2005) simulated the Storegga slide as a series of retrogressive blocks on a 2.08 Â 2.08 km grid for the Norwegian-Greenland sea, with a nested 500 Â 500 m grid focused on a limited region of the Norwegian coast. This work was extended by Løvholt et al. (2005) to include ideas about how the slide may have moved. A major limitation of these studies was an inability to resolve complex coastlines in the regional models, hence the use of nested models. In particular, no study to date has quantified the effect of increasing coastline resolution on the numerical simulations. An alternative to nested models is to use a multiscale simulation, where grid resolution varies spatially, often by orders of magnitude . Multiscale models often use an unstructured mesh, so in addition can accurately represent complex coastlines and bathymetry without ''staircase'' effects (Wells et al., 2005). Multiscale modelling then also allows more complex coastal morphologies to be included in the simulation.
Here, we use Fluidity-a 3D finite element, non-hydrostatic, numerical model that makes use of unstructured triangular/tetrahedral meshes to enable accurate representations of the domain and allow multiscale simulations of large regions. Fluidity has previously been used to simulate earthquake-generated tsunami (Shaw et al., 2008;Mitchell et al., 2010;Oishi et al., 2013). Oishi et al. (2013) showed that Fluidity could accurately simulate the 2011 Japanese tsunami and, in particular, was able to represent the dispersive effects of the tsunami by using multiple vertical layers. We do not consider the effects of dispersion here due to the size of the Storegga slide. The advantage of our non-hydrostatic model over widely used non-dispersive shallow water models is that it can be used to capture wave dispersion by including multiple vertical layers (Oishi et al., 2013, e.g.), but can also approximate the shallow water approach using a single layer to model the propagation of non-dispersive waves (Mitchell et al., 2010, e.g.). As the slide-tsunami scenarios we investigate here generate non-dispersive or very weakly dispersive waves our simulations generally use only a single layer. While this results in a (modest) computational overhead compared to alternative formulations, the benefit is that the results presented here can be directly compared with future studies, using the same model, that examine highly dispersive waves generated by, for example, smaller slides (Glimsdal et al., 2013, e.g.). Mitchell et al. (2010) used the same model to study ancient tsunamis in the Jurassic Tethys sea, which shows the flexibility of the model in representing arbitrary coastlines. Here we describe how Fluidity has been modified to simulate slide-tsunami generation using prescribed rigid-block slide motion. This allows two of the four phases of slide-generated tsunami waves to be studied : the generation and propagation of the wave. The simulation of slide dynamics and tsunami wave inundation are not considered in this work.
Previous studies of the Storegga slide tsunami did not directly include the effects of relative sea-level changes on bathymetry (Harbitz, 1992;Bondevik et al., 2005). Isostatic adjustments from ice-sheet loading and unloading produce complex changes in relative sea-level across the region. Recent studies have simulated this process to produce 1000-year time slices of such changes since the Last Glacial Maximum (Bradley et al., 2011). Relative sea-level changes of up to 50 m have occurred since the Storegga slide, which caused substantial changes in coastlines. For example, 8000 years ago a region in the southern North Sea was an island-Doggerland (Fitch et al., 2005)-and the coastlines around Norfolk, UK, and the northern coast of mainland Europe (Fig. 1) were dramatically different. Human artefacts (flints and spearheads) and mammal remains (mammoth and rhinoceros teeth) have been dredged from the Dogger Bank (Flemming, 2002). There has been speculation that the Storegga tsunami was the cause of the abandonment of the island by Mesolithic tribes (Weninger et al., 2008).

Aims
In this paper, we first briefly describe the Fluidity model and the newly-implemented rigid-block slide model used to initiate the tsunami. We verify the implementation of this model by comparing our results to previous numerical results for test problems in both 2-and 3-dimensions. We then demonstrate numerical convergence of the model, before using a multiscale model to show how it is possible to achieve high resolution in coastal areas in a simulation that includes the entire Norway-Greenland sea and tracks the tsunami wave for 15 h. Our aims are to demonstrate the effectiveness of multiscale simulations for slide generated tsunamis. Finally, we show the effect of incorporating palaeobathymetric changes on the simulated run-up heights.

Fluidity
Fluidity is a highly flexible finite-element/control-volume modelling framework which allows for the numerical solution of a number of equation sets  and has been used in a number of flow studies ranging from laboratory-to oceanscale (e.g. Wells et al., 2010;Hill et al., 2012;Hiester et al., 2011). In an ocean modelling context, Fluidity has been used to model both modern and ancient earthquake-generated tsunamis (Oishi et al., 2013;Mitchell et al., 2010;Shaw et al., 2008). Here, Fluidity is used to solve the non-hydrostatic incompressible Navier-Stokes equations under the Boussinesq approximation in a rotating reference frame: @u @t where u is the 3D velocity vector, t represents time, p is pressure, m is the kinematic viscosity tensor (isotropic and set to 1 m 2 s À1 ) and q denotes the density, which is constant in this work. The reason for choosing an isotropic viscosity is that experiments showed no discernible differences in results when using different values of viscosity in the horizontal and vertical when using a single layer of elements. This may not be the case when multiple layers are used to capture dispersion (Oishi et al., 2013). X is the rotational velocity of the Earth and g is the gravitational acceleration with k pointing in the radial, upward direction. Eq. (1a) is discretised using a linear discontinuous Galerkin approximation (P1 DG ) for velocity. A pressure projection method is used to solve for the pressure p and enforce a divergence-free velocity field at the end of each time-step. Pressure is discretised using a continuous Galerkin, piecewise quadratic formulation (P2). The resulting P1 DG P2 velocity/pressure discretisation has a number of desirable properties described fully in Cotter et al. (2009a,b) and Cotter and Ham (2011). A two-level h method is employed for time-integration. Here h ¼ 0:5 which yields a second-order accurate, implicit Crank-Nicolson scheme. Two Picard iterations per time-step are used to linearise the nonlinear advection term.
A combined pressure-free-surface kinematic boundary condition formulation is employed as the top boundary condition (Funke et al., 2011;Oishi et al., 2013). A no-normal flow with a quadratic bottom drag, with dimensionless coefficient C D set to 0.0025, is applied at the bottom, except where the slide motion is prescribed (see Section 2.2). At the coastlines a free-slip no-normal flow formulation is used and at the open boundaries either a velocity or a free surface elevation is prescribed. Further details of the discretisation methods employed are given in Piggott et al. (2008) and AMCG, Imperial College London (2014).

Slide motion
The Storegga slide was a large submarine slide which disintegrated during movement (Haflidason et al., 2005), such that it was not a single rigid block. Moreover, there is evidence that slope failure started in deep water and moved retrogressively upslope (Masson et al., 2010). However, as such complex slide dynamics would add considerable computational expense, here we adopt a simplified slide movement formulation described by Harbitz (1992) and Løvholt et al. (2005). The slide is a rigid block that has a prescribed shape and moves using a prescribed velocity function. Despite its simplicity, Storegga-tsunami simulations using this approach produced run-up height estimates in reasonable agreement with those inferred from sediment deposits at a range of locations (Bondevik et al., 2005).
The total water displacement is determined by the changes in aggregated thickness as the slide moves with a prescribed velocity. We impose this water displacement as a normal velocity Dirichlet boundary condition, u Á n ð Þ D , calculated as: u Á n ð Þ D ¼ À h s ðx À x s ðt À DtÞ;y À y s ðt À DtÞÞ ½ À h s ðx À x s ðtÞ; y À y s ðtÞÞ ½ Dt ð2Þ where Dt is the timestep of the model, and n is the outward unit normal. The slide motion is defined as: hðx; y; tÞ ¼ h s ðx À x s ðtÞ; y À y s ðtÞÞ; where hðx; y; tÞ is the slide thickness in two-dimensional Cartesian space ðx; yÞ at time, t, and h s is the vertical displacement (with respect to the boundary) of water by the slide.
The parameters x s and y s describe the slide motion and h s describes the slide shape via simple geometric relationships: Here, / is the angle from the x-axis that the slide travels in, ðx 0 ; y 0 Þ is the initial position of the centre of the slide front, R is the run-out distance, and, T is the total time of the slide travel, defined as: where T a is the acceleration phase of the slide, T c is the constant speed phase, and T d is the deceleration phase. The acceleration time T a ¼ pR a =2U m (acceleration distance R a ), the constant speed time , and the deceleration time , define the relationship between travel time, maximum speed, and run-out distance for the three phases. The total run-out distance of the slide is , where an overlay of the production mesh used in this study is also shown. Shading shows water depth with darker shades indicating deeper water. For the insert the modern coastline is also shown (light grey) over the palaeo-coastline (dark grey).
The term sðtÞ in (4) governs the acceleration and deceleration phases, given a maximum slide velocity U max , and is defined as Acceleration phase: Constant speed phase: Deceleration phase: The slide shape is defined as: where the slide has dimensions of maximum height, h max , length, L, and width, B. To avoid sharp edges, which would cause numerical oscillations, a smoothing length, S, is used at the front and back of the slide and the slide is smoothed along the whole width laterally as described in Harbitz (1992). S is 1 km in the 2-D validation study and 7.5 km in the Storegga simulations. The slide movement is then governed by x 0 and y 0 , which describe the slide motion in the x-y plane and are defined by: This gives a total volume of the slide, V: The motion given by (2) is then weakly imposed in the normal direction on the lower boundary to simulate the rigid block slide. This is a similar method to Ma et al. (2012) and Harbitz (1992), though differs in that Harbitz (1992) alter the h term in the shallow water equations. In practice, all methods should give very similar results.

Model validation
To ensure correct operation of the slide-tsunami model for weakly dispersive or non-dispersive waves we replicated simulations from independent numerical modelling studies in the literature. The first is a flat two-dimensional model, with dimensions approximately equivalent to the Storegga slide , which produces a non-dispersive wave. The second is a smaller-scale, three-dimensional slide on a gentle slope (Ma et al., 2013), which produces a weakly dispersive wave. Comparisons to these previous studies verify correct implementation of the slide boundary condition. Fluidity's ability to capture highly dispersive slide-tsunami will be examined in future work. Haugen et al. (2005) simulated wave generation by the Storegga slide using a two-dimensional (x-z) approach, with an idealised rigid-block slide geometry and constant water depth. They showed that the very large length of the Storegga slide compared to the water depth resulted in a very long wave with little-to-no dispersive characteristics. Here we reproduce this simulation using Fluidity with a single element in the vertical. Tests with more vertical layers (not shown) produced almost identical results, confirming that wave dispersion is negligible in this scenario. The test case uses a flat-bottom domain, 1000 m deep, and 2000 km long. The slide has the parameters detailed in Table 1. Fluidity simulated the same scenario at six different horizontal resolutions: 5000 m, 2000 m, 1000 m, 500 m, 250 m, and 125 m. The mesh in this case is formed of 1D elements in the horizontal, which are then extruded downwards to 1000 m. A single layer of triangular elements was used in the vertical and the timestep was fixed at 1 s. The Fluidity results are compared to Haugen et al. (2005) in Fig. 2. Qualitatively, the results are identical at high resolution, with consistent peak amplitudes that occur at the same locations. The numerical oscillations visible in the Fluidity output become negligible at 1000 m resolution, with little difference between results at resolutions between 1000 m and 125 m (Fig. 2). The observed numerical oscillations are caused by the sharpness of the leading and trailing edges of the slide, where minimal smoothing of 1000 m was used . Increasing the smoothness of these edges (by increasing S in (9)) removes the oscillations. Clearly, the mesh resolution must be high enough to capture the smoothing length or the slide will have an effective flat front. To check that this was the cause of the spurious oscillations, the 5000 m resolution case was re-run with a smoothing length of 7500 m. The results show much reduced oscillations, but with the wave form shifted due to the new location of maximum height (Fig. 2). This experiment confirms the correct implementation of the boundary condition and shows how the assumed shape of the slide dictates the mesh resolution required in the slide area. A slide with steeper leading and trailing edges requires higher spatial resolution to eliminate numerical oscillations.

3D validation test case
To extend our validation of Fluidity's new slide-tsunami model to three dimensions, we also replicated a simulation of landslide generated waves that are only weakly dispersive (Ma et al., 2013). Recent work by Ma et al. (2013) simulated the wave train produced by a rigid-block model in a three-dimensional domain on a constant slope. We can therefore compare Fluidity to the results shown in Ma et al. (2013). The domain is 8 Â 8 km, with a constant slope of 4°. We set the minimum depth to be 12 m and the maximum to be 400 m. We used a horizontal model resolution was 25 m in x and y and explored the influence of vertical resolution by performing simulations with 1-4 layers. Ma et al. (2013) use a different slide geometry to that described above, based on the work of Enet and Grilli (2007). The slide geometry is given by: where The truncation parameter, is 0.717.
The slides moves according to: where s 0 ¼ u 2 t =a 0 ; t 0 ¼ ut a 0 ; a 0 ¼ 0:27 m s À2 , and u t ¼ 21:09 m s À1 as detailed in Ma et al. (2013). We use these definitions of the slide height and speed for comparisons to Ma et al. (2013).
The resulting wave is very similar in magnitude and waveform to that shown in Ma et al. (2013), even using only a single layer in the vertical (Fig. 3). Convergence of the Fluidity model results is observed for three or more element layers (c.f. 40 layers used by Ma et al. (2013)), indicating that the wave is only weakly dispersive. In more detail, Fluidity produces slightly lower amplitude waves than those reported by Ma et al. (2013) (Fig. 3), at earlier times, although Fluidity then produces higher positive amplitudes by 100 s. It is not clear if the model described by Ma et al. (2013) overestimates wave height or Fluidity underestimates. It should be noted that previous comparisons of Fluidity to both numerical models and observational data, Haugen et al. (2005) and Oishi et al. (2013), show excellent agreement to both amplitude and phase of wave patterns resulting from both slides and earthquakes in two-and three-dimensions at ocean scales.

Model set-up
Having benchmarked the implementation of the prescribed slide boundary conditions against independent models, we now show how Fluidity is capable of simulating real-world scale slide-generated tsunamis with high resolution in areas of interest by recreating the Storegga slide.
The same domain is used for all simulations described here. The domain stretches from 43°west to 24°east and 47°north to 80°n orth. GSHHS data (Wessel and Smith, 1996) was used to generate coastlines for all modern simulations, which has resolutions of 200 m (full) to 25 km (coarse). For the simulation involving palaeobathymetry the coastline was derived from the 0 m contour. Bathymetric data was derived from GEBCO (IOC, 2008) which has resolution of 1 arcminute (approximately 2 km in this region). For each domain QGIS (QGIS Development Team, 2009) was used with bespoke software to generate coastline input for GMSH  Haugen et al. (2005). As the Fluidity mesh is refined, numerical oscillation cease to be observed and the result become very similar from that of Haugen et al. (2005). Numbers given in the key are the mesh resolution in metres. (b) Free surface height variation at time 114 s after slide initiation for an idealised Storeggatype slide, based on the numerical experiment of Haugen et al. (2005). Smoothing the mesh reduces numerical oscillations as expected but does slightly shift the timing of the wave. (Geuzaine and Remacle, 2009) which created the horizontal computational mesh. The mesh is on a Cartesian sphere of radius 6371.01 km. Coastlines were constructed using a B-spline curve through the points given by the GSHHS data. Bathymetry is incorporated by extruding the generated surface mesh radially downward to the depth given by the bathymetric data, which is carried out at run-time. Each simulation uses a one-element deep solution, effectively a depth-averaged velocity as used in (Mitchell et al., 2010;Wells et al., 2010). A consequence of this approximation is that a minimum water depth has to be specified for the mesh as inundation (wetting and drying) was not utilised in this study. Here, a minimum depth of 10 m was used.
We generate the slide using the single rigid block slide, described in Eqs. (4)-(11), following the work in Harbitz (1992), using the parameters in Table 2. Note that we do not include the effects of retrogressive slide evolution. This style of multi-block slide motion was investigated in Løvholt et al. (2005) and Bondevik et al. (2005), who concluded that the time interval between block initiation would need to be very small in order to produce large wave heights consistent with observation and such scenarios are qualitatively similar to the motion of a single continuous body.
For initial runs, to explore the sensitivity of model results to spatial resolution, the simulation was run for five hours model time, which was sufficient to allow comparison with previous studies. A preliminary study of model sensitivity to time step choice found little difference between a time step size of 1 s and 2 s, and hence 2 s is used in all simulations. Even though the Crank-Nicolson discretisation provides stability at larger time steps, with implicit schemes one still needs to select an appropriate time step size to ensure model accuracy. Here we are concerned with the propagation of waves over relatively large distances and the implicit discretisation employed here tends to damp these waves if too large a time step size is used, see also Oishi et al. (2013). The robustness that comes with the use of implicit time stepping schemes is particularly useful when an unstructured mesh of a complex region might include particularly small elements in order to resolve complex coastlines, and these could significantly impact on the time step restrictions with a fully explicit model. A similar issue arises in the use of flooding models (which will be considered in future work) where the inundation front may propagate large horizontal distances in very short time scales (Funke et al., 2011). The final two simulations using multiscale resolution were run for 15 h to track the wave propagation as far as Doggerland and the English Channel.  Previous studies of the Storegga slide-tsunami have not included the changes in bathymetry that have occurred in the last 8.15 kyr (Harbitz, 1992;Bondevik et al., 2005). We test the effect of that in this work. In addition, model predictions of wave heights are also sensitive to slide geometry, retrogressive behaviour, acceleration and maximum speed. These will be explored in future work; first we have to establish confidence in the numerical factors. We compare results to the virtual wave gauge records shown in Harbitz (1992) and Bondevik et al. (2005), which in turn are compared to inferred run-up data, where available. The location of these gauges are shown in Fig. 4. Harbitz (1992) used eight wave gauges placed around the Norway-Greenland sea and in the vicinity of the Storegga slide. Bondevik et al. (2005) detail 25 sites where run-up heights can be estimated and show the free-surface variation with time at seven of these locations. We added a further two gauges on the east coast of Scotland (26) and north-east England (27). In addition, Bondevik et al. (2005) performed an experiment where they varied resolution in a small subdomain around Sula, Norway and showed the effect of resolution on simulated wave height observed there. We compare Fluidity against the highest resolution (500 m) results given by Bondevik et al. (2005).

Mesh resolution
To examine the effects of horizontal mesh resolution the domain was constructed using the coarsest resolution GSHHS coastline data, which has a resolution of around 25 km. A constant element edge length was then defined to match 50 km, 25 km, 12.5 km, and 6.25 km. No mesh metric was used to alter mesh based on, for example, distance to coastline, and hence the meshes had the same resolution across the whole domain. Note that the resolution of the coastline was such that all mesh resolutions were high enough to resolve almost all coastal features. For the 50 kmmesh, a few small islands could not be properly resolved due to the higher resolution of the GSHHS data and these islands were therefore removed from the meshes at all resolutions. The number of tetrahedral elements changed by a factor of approximately four for a doubling in resolution, such that the 6.25 km resolution simulation contained nearly 64 times the number of elements as the 50 km resolution simulation (Table 3). Due to the increase in element count, the modern multiscale simulation was carried out on 540 cores on the Imperial HPC system. Run time was approximately 56 h for 15 h of simulated time. Note that no parallel scaling tests were performed to ensure maximum parallel efficiency. In general runtimes are proportional to the number of elements, which in turn in proportional to the number of degrees of freedom.
In addition to discretisation errors, the change in resolution has two consequences. One is an improvement in the representation of bathymetric data by the computational mesh and the second is a change in the position of the virtual wave gauges (see Section 4). The bathymetry used here is the GEBCO 1 arcminute data. This is equivalent to $1.8 km resolution at this latitude, so even the highest resolution mesh used in this mesh resolution experiment cannot resolve all bathymetric features. We interpolate the bathymetry to each vertex in our computational mesh using bi-linear interpolation. As the mesh is refined, more features are resolved. The second effect is the refinement of detector locations. In order for a detector to be contained with the mesh (i.e. not on land as represented by this coastline), the latitude and longitude position was converted to spherical Cartesian coordinates and the detector was then moved to the closest mesh vertex. Similar issues occur in other studies (Bondevik et al., 2005).

Multiscale run
The multiscale mesh (Fig. 5) was constructed in a similar manner to those above. Resolution varied from 500 m to 50 km and resolution was dependent on bathymetry, Hessian (second-order gradient matrix) of the bathymetry, distance to coastline (see Lambrechts et al., 2008 for details) and distance from slide location. Distance from slide was determined by tracing the approximate slide locations through time and then using GDAL (GDAL Fig. 4. Location of virtual wave gauges taken from Harbitz (1992) and Bondevik et al. (2005) and the two additional gauges used in the simulations presented here. Shading shows the GEBCO bathymetry and grey solid region shows land areas of the modern multiscale simulation. Development Team, 2013) to generate a mesh with resolution of 2 km in a region in the slide area, and which smoothly increased to a mesh spacing of 50 km at 100 km distance from the slide region. Coastlines were generated from GSHHS. The UK, Ireland and neighbouring islands were generated using the full resolution dataset (which has an approximate 200 m resolution). All other coastlines were generated using the intermediate resolution data (which has an approximate 1 km resolution). Small, unresolved islands were removed from all coastlines. Due to the different coastline resolutions, the UK and Irish coastlines were meshed using 500 m element edge lengths, whereas 1 km element edge lengths were used along other coastlines. Note the maximum mesh element size was very coarse, but this only occurred in very deep, flat areas away from the slide region: all shallow regions or regions where depth varies rapidly had much finer resolution due to the choice of metric. Note also that the horizontal resolution around the coastlines was much less than that of the bathymetry data (1 km or 500 m mesh resolution vs. 1.8 km bathymetry resolution) and hence all bathymetric features were well resolved in these regions.

Palaeobathymetry run
The palaeobathymetric domain was generated by first adding the isostatic adjustment data from Bradley et al. (2011) to the GEBCO bathymetry dataset to generate a palaeobathymetry. Note that the isostatic data only has extent of À20°to 20°west to east and 40°to 70°south to north, and hence we extrapolated the data by setting the extended domain corners to the same values as the corners of the true domain and then using GMT to interpolate the missing data. We extrapolated data to match our domain (À43°to 24°west to east and 22°to 80°south to north). Results from this simulation are therefore only valid within 20°west to 20°east and 43°north to 70°north. Note that all wave gauges are situated within this region except gauge 1 (Greenland). All comparisons to the multiscale mesh were carried out within this sub-domain. Once the palaeobathymetry was generated, the 0 m contour was used to generate a coastline as GSHHS was no longer valid. Inland seas and lakes were removed. Mesh resolution, including refinement in the vicinity of the slide and around bathymetric features, was identical to the modern multiscale simulation, except all coastlines were generated using 1 km element lengths. As before any small islands and features were removed if they could not be resolved. The resulting coastline and bathymetry are shown in Fig. 1 which also shows the comparison to the high resolution GSHHS data. There are clear differences in coastline configuration around the eastern coast of the UK, but no significant differences around the central and southern Norwegian coasts. The mesh contains just over 1 million elements, around 300,000 fewer than the modern mesh, which is largely due to the difference in coastline resolution and the reduced ocean area (Table 3).

Results
For each simulation we compare the basin-wide free-surface (i.e. sea surface) height and the free-surface variation at the 34 virtual wave gauges. We compare against a subset of these locations for each simulation. Fig. 6 shows the large-scale free-surface patterns and the qualitative convergence between 25 and 12.5 km mesh resolution. There are no discernible differences in free surface at 60 min simulated time for resolution of 25 km and below. Minor differences between the 25 km and 12.5 km simulation output at 120 min can be seen, but there is no visible difference between 12.5 km and 6.25 km resolution output at both 60 and 120 min (Fig. 6). The 50 km resolution is clearly too coarse to provide accurate information on the wave form, but is adequate to provide information on first arrival times. Note that the slide smoothing parameter, S is 75 km, which is less than two mesh elements in size for the 50 km resolution simulation. This is likely to be the primary cause of the numerical oscillation observed.

Mesh resolution
Examining the results in more detail shows some differences between results from simulations with 12.5 km and 6.25 km mesh resolutions (Fig. 7). Using 50 km mesh resolution often leads to numerical oscillations in the solution, with peak wave heights that are out-of-phase of the higher resolution simulations. These resolutions are caused by the smooth slide edges not being resolved correctly. Similar oscillations can be seen at 25 km mesh resolution at some locations (e.g. gauge 4) and can also show anomalously large wave heights, for example at gauge 9. Once mesh spacing is below 12.5 km, these oscillations do not occur, and for many locations the difference between 12.5 km mesh resolution and 6.25 km mesh resolution is relatively small. We can therefore conclude that 12.5 km mesh resolution is suitable to minimise numerical effects on the solution. In addition, this also gives a reasonable number of elements in the computational mesh (Table 3).

Multiscale modelling
From the experiments described it is clear that the large-scale simulated results do not depend on bathymetric data sources or mesh resolution once numerical convergence has been achieved. However, it is also clear that at coastal-scales the resolution of the bathymetry and coastline can alter the results obtained considerably, often in non-intuitive ways. An obvious solution to this issue is to use multiscale resolution where the resolution across the majority of the domain can be low and then be refined over areas of interest, coastlines and around changes in bathymetry. Using the multiscale mesh described in Section 4.2, we performed a 15 h simulation of the Storegga tsunami using an otherwise identical set-up to that described in Section 2.2. We compare the results to estimated run-up measurements from observations as well as previous results above. Note that the mesh is large, containing some 1,378,146 elements (Table 3), which is around 300,000 fewer than the fixed mesh 6.25 km resolution simulation. The number of elements can be reduced further by reducing the coastline resolution around the UK and around the Storegga slide itself which should result in little difference to the results presented here. Further work is required to optimise the mesh for computational efficiency without loss of accuracy.
Results from this simulation are similar to those in previous experiments in the observed free-surface variation at both one and two hours. There is also an expected reduction in maximum Fig. 6. Simulation of the tsunami at 120 min after slide initiation. The different resolutions are 50 km (a), 25 km (b), 12.5 km (c), and 6.25 km (d). There is very little difference in the output once resolution reaches 25 km.
wave height at several locations as a result of increased bathymetric and coastline resolution, along with an increase at other locations. All observed features in all wave gauges are consistent with the behaviour seen in the above experiments. The basin-scale free-surface variations are indistinguishable from the 6.25 km resolution simulation (Fig. 8).

Comparison to deposits
Observational run-up height estimates of the incident wave of ancient tsunamis are inferred from the location of high-energy sedimentary deposits that can be traced inland or between raised lakes (e.g. Bondevik et al., 2005). Such estimates are generally underestimated as this is the minimum run-up height required to explain the deposits. For the Storegga slide there are a number of observations in northern Scotland and along the Norwegian coast, as well as one mapped deposit on the Faroe Islands. The maximum simulated wave height can be compared to inferred wave heights at these locations. Fig. 9 shows the free-surface heights at key locations where tsunami deposits have been found, with estimates of the run up heights included following Bondevik et al. (2005). Note that the fixed horizontal resolution of 6.25 km does not always match the multiscale resolution results, e.g. gauges 24 and 12 (Fig. 9), highlighting the need for high resolution in coastal regions . For the multiscale simulation there is good agreement at all stations, with exception of those around the Faroe Islands (32) where our models (and those of Bondevik et al., 2005) underestimate the wave height. A good agreement with estimated wave heights is found at Sula, Norway (15), where Bondevik et al. (2005) simulated a 20 m wave, but estimated a 10-12 m from sediment deposits. Our models predict a wave height of 14.5 m, which is a better agreement. Similarly, Brønnøysund and Hommelstø in northern Norway (wave gauges 9-11) have an estimated run-up height of >3 m (Bondevik et al., 2005), but previous simulations predict a 17.9 m wave, which is probably a large overestimation (Bondevik et al., 2005). Here, we  Harbitz (1992) and Bondevik et al. (2005) to fluidity at resolutions from 50 km to 6.25 km. record a maximum wave height of 5.8 m, which is a more reasonable result. Around the Shetlands we predict a wave height of around 8 m, lower than that estimated from deposits, but an improvement on previous modelling efforts (Bondevik et al., 2005).
The results using palaeobathymetry show little difference to those using modern bathymetry except at a few key locations. The large-scale features north of 52°N show very little difference (Fig. 8). The maximum wave height in the domain is largely unaf-fected by the inclusion of palaeobathymetry (Fig. 10), with most of the study area experiencing a difference in wave heights of only a few metres. However, smaller regions show a substantial increase in maximum wave heights, in particular the Shetland Islands, where maximum wave height increases by nearly 5 m when using palaeobathymetry (Fig. 10g). This gives an improved match to estimated run-up heights, which were several metres too low in previous studies (Bondevik et al., 2005). An area off the northern Fig. 8. Snapshots of free surface 2 h after slide initiation. Left (a) uses the coarse GSHHS coastline data with constant 6.25 km resolution, right (b) shows the results from the multiscale simulation which includes high and full resolution GSHHS coastline data. Bottom (c) simulation using palaeobathymetry. There is little difference at this scale to the simulated wave form. Note that land mask shows the full detail of the simulated coastline in each simulation result. Fig. 9. Free surface variations at key virtual wave gauge localities where run-up heights can be estimated from sediment deposits. These are concentrated on the coasts of Norway and Scotland. Plots of free surface variation from both the multiscale and palaeobathymetric simulations are shown. Where available, the equivalent virtual wave gauge from Harbitz (1992) or Bondevik et al. (2005) is also shown. Numbers refer to detector locations shown in Fig. 4 coast of Norway shows a similar increase in maximum wave heights, though the difference is small once the wave reaches the coastline (Fig. 10c). The north-eastern coast of the UK experienced waves between 3-6 m, much like the eastern coast of Scotland, although only one possible deposit has so far been found (Boomer et al., 2007). The southern North Sea, especially the coasts of the UK and Dogger Bank show significant differences, largely due to the alteration of the coastline, but there are no known observations here. Wave heights are predicted to be around 1 m on the UK coast and up to 5 m on the northern coast of Doggerland. The maximum elevation of Doggerland here is less than 10 m, with large areas of less than 5 m. It is therefore possible that much of Doggerland would have been flooded by such a wave. Due to the inclusion of the Doggerland island, the northern coast of mainland Europe experiences maximum wave heights of 1 m or less -much lower than if modern bathymetry is used. The wave also reaches the western coast of the UK, with maximum wave heights of around 1 m on the Cornwall and Devon coasts. Similarly we predict waves of up to 5 m on the western coast of the Republic of Ireland.
On a more local scale locations such as gauge 7 show a significant shift in the arrival time of the waves (9). Many locations show a slight increase (e.g. 30) of a few metres, which improves the match to estimated run-up heights (9), whilst a number show very little difference (e.g. 15). All other locations where Storegga tsunami deposits are found show a good match to observed data using either palaeo-or modern bathymetry, with the exception of the Faroe Islands where the wave height is underestimated and the inclusion of palaeobathymetry makes little difference. The modern result is very similar to that of Bondevik et al. (2005) who postulate that the wave is amplified in the fjord. We therefore conclude that palaeobathymetry can have a significant effect at a local scale, similar to the increase in bathymetric and coastal resolution, but has little effect on the basin-scale results. We also note that at some locations, such as the Faroe Islands there is little difference in the modelled wave height, despite a significant drop in relative sea level of around 20 m in the region. However, the changes in relative sea level also affect the propagation of the wave along the wave path to the Faroe Islands, so it is overly simplistic to use the modern bathymetry and account for the change in relative sea level at a single location. The discrepancy here may be due to local funnelling or amplification effects and a further increase of resolution may resolve this.
Videos of these two simulations are available in the supplementary material.

Multiscale modelling
The idea behind multiscale resolution simulations is that areas of interest can be simulated at an appropriate resolution without the expense of computational effort in areas where high resolution is not required. Nested models can also be used to similar effect. Allgeyer et al. (2013) used a series of nested finite-difference grids to examine the effect of the Lisbon 1755 tsunami on tidal gauges in La Rochelle, France. Grids were nested from 1 0 (approx. 2 km) to 0.3 00 (9 m), zooming in on the target region. No sensitivity to mesh resolution was carried out, however. In addition, Roger et al. (2010) used the same method to study the effect of the Lisbon 1755 tsunami on Caribbean Guadeloupe Archipelago, with similar resolutions to Allgeyer et al. (2013). These two studies nested the same computational model; however, it is also possible to nest different models to carry out large-scale simulations. Kirby et al. (2013) used the non-hydrostatic model of Ma et al. (2012) in the nearfield source domain, before linking this to the larger-scale model described in Kirby et al. (2013) to investigate the 2011 Japanese tsunami. Resolution varied from 1 km to 2 0 (approximately 4 km). Horsburgh et al. (2008) followed a similar methodology to study the effect of the Lisbon 1755 tsunami on the UK coast, using a finite-difference model with approximately 3.5 km resolution in the larger domain and a finite-element model around the UK coast with resolution varying from 10 down to 1 km. It is clear with all of these studies that resolution around areas of interest is important, Fig. 10. Maximum sea surface height (positive wave height) and difference maps. Difference maps are modern and palaeobathymetry maximum free surface height. The maximum wave height for the simulation using modern bathymetry (a) and palaeobathymetry (b) show little difference. However, plotting the difference between these (c) shows difference of up to 5 m around Norway, the southern North Sea ((d) -maximum sea surface and (e) -difference map), and the Shetlands ((f) -maximum sea surface and (g) -difference map). For the difference maps negative numbers indicate where the maximum free surface is less using modern bathymetry (i.e. palaeobathymetry wave heights are greater) and positive numbers indicates the maximum sea surface is greater using modern bathymetry (i.e. palaeobathymetry wave heights are lower). but all must limit their regions of interest. The multiscale modelling technology shown here can allow multiple areas of interest within the same simulation, whilst capturing changes in bathymetry and coastline in the mesh. It is also worth noting the lack of studies detailing the effect of resolution for tsunami simulations. Bondevik et al. (2005) did show a clear convergence of results using a smaller region simulation at both 250 and 500 m resolution. The technology presented here could be further improved by increasing resolution even further to that used by other studies above, for example 10 m, around a particular small region of interest.
As part of this work we investigated the effect of a number of factors on the estimated run-up heights of the tsunami. These were: bathymetric data source (GEBCO or ETOPO (Amante and Eakins, 2009)), the resolution used to generate the coastlines and the bathymetric resolution. From these experiments only coastline resolution made a substantial difference. Virtual wave gauge 24 (Fig. 9) shows an example where the effect of coastline resolution makes a substantial difference to the estimate run-up height as the high resolution fixed mesh case (using the coarse resolution GSHHS data) produces a much large wave height than the multiscale mesh where the high resolution GSHHS data were used. There are also virtual wave gauges (not shown) that show an increase in wave height with increasing coastal resolution. However, it is difficult to ascertain if it is increase in coastal resolution or the associated increase in bathymetric resolution due to the increased mesh resolution that is the primary cause of these changes.

Model limitations
There are some limitations in the present study. The lack of inundation at the coastlines, coupled with the minimum depth requirement, means that the true free-surface variation at an arbitrary coastal location cannot yet be represented. Fluidity is capable of simulating inundation in a limited region (Funke et al., 2011) and work is ongoing to link this technology to large-scale simulations. The virtual wave gauges must be contained within the mesh to record the free surface variations at a given location. As we varied coastlines and resolution, wave gauges were moved slightly between simulations to ensure they were not on land. Bondevik et al. (2005) used a similar methodology as the gauges specified there were not within their computational domain. They do not report the true location as the effect of this shift was thought to be small. The largest difference in the present study was less than 1 degree for the 50 km resolution simulation with the coarsest GSSHS coastline. All other simulations had differences of much less than 1 degree.
The current model does not include inundation as the wave reaches the coastline. Therefore comparisons are made between the estimated run-up height from sedimentary deposits and the maximum wave height in the vicinity of the deposit. The difference between the two estimates will depend on local factors, such as vegetation and small-scale (i.e. unresolved) bathymetric/topographic changes. We aim to include this in future work.
Perhaps the most important simplifying assumption within this study is that the Storegga Slide moved as a single rigid block. This a priori assumption is important because the way in which the original slide moves determines the initial dimensions of the resulting tsunami. Field observations (Haflidason et al., 2005) suggest that much of the slide mass disintegrated, such that it was not a single rigid block. Moreover, there is evidence that slope failure started in deep water and moved retrogressively upslope (Masson et al., 2010). This modelling also assumes a priori that the slide accelerated to a speed of $35 m/s over 3365 s. The acceleration trajectory of the slide is unknown, although previous modelling suggests that such fast speeds are needed to generate a large far field tsunami.
We have based our model on the work of (Harbitz, 1992). This was later refined in terms of both the slide shape and initiation by Bondevik et al. (2005) but no comparison to Harbitz (1992) was carried out and hence it is difficult to ascertain what effect these modifications had on the model results. Bondevik et al. (2005) do not give an analytical expression for the modified slide and hence it could not be used in this study. In addition, Bondevik et al. (2005) also increased resolution of the mesh from 12.5 km to 2.08 km, possibly confounding any comparison. This is the first study of the Storegga tsunami, to report on the effects of mesh, coastline and bathymetric resolution on a simulated slide-initiated tsunami. We show a good match to estimated wave heights, but these might be further refined by adjusting the slide parameters further, as per Bondevik et al. (2005). The Fluidity modelling presented here assumes one particular type of slide movement as a single rigid block. It is unclear how somewhat more realistic slide behaviour would affect tsunami magnitudes and inundation heights around surrounding coastlines. More work is required in order to attempt to improve the veracity of the model by altering the slide initiation and shape and to study the effects of such changes and how they compare to the changes described here with respect to resolution.

Conclusions
The effects of bathymetric and coastline resolution are important in determining accurate simulated run-up heights of tsunamis. We have shown that the higher resolution coastline and bathymetric simulations produce simulated wave heights that are in closer agreement to inferred wave heights from observational data and have some sense of numerical convergence. Overall numerical resolution is important to minimise numerical errors and for this simulation a fixed mesh of 12.5 km is sufficient with coarse coastlines to reproduce the work of Harbitz (1992). However, as along-coastline resolution increases, commensurately higher mesh resolution is required around the coasts.
Assumptions of the slide acting as a rigid block, accelerating to 35 m/s, are similar to previous studies, but as the Storegga slide is thought to be retrogressive and disintegrate as it moved, more work is required to ascertain the effects of this on wave run-up heights. In establishing the spatial resolution of coastlines and palaeobathymetry required to adequately model the Storegga slide-generated tsunami, this work provides a foundation on which simulations examining the effect of complex slide parameters can build.
Given the simplicity of our slide model and the absence of an inundation model, our multiscale models of the Storegga submarine slide generated tsunami shows remarkable agreement with inferred wave-heights from sediment deposits along the Norwegian and Scottish coasts. The agreement within the Faroe Islands is less good, with a simulated wave height that is around a 6 m too small, but consistent with previous studies (Bondevik et al., 2005). Our multiscale model simulates the Storegga tsunami for 15 h, tracking the wave propagation into the southern North Sea, predicting wave heights of less than 1 m for the northern coast of mainland Europe. The addition of palaeobathymetric information, neglected in previous studies, aids the match to observed data within the region where our data is valid and makes a substantial difference in the southern North Sea region and around the Shetland Islands. However, the use of realistic palaeobathymetry makes little difference along the Norwegian coast, which was the primary focus of previous studies. As an example of the importance of considering palaeobathymetry, we show that Doggerland would have experienced wave heights of up to 5 m. Given the majority of this island was less than 5 m in height, it would have experienced wide-scale flooding. It is therefore plausible that the Storegga slide was indeed the cause of the abandonment of Doggerland in the Mesolithic.