Reconciling Storegga tsunami sedimentation patterns with modelled wave heights: A discussion from the Shetland Isles field laboratory

The Shetland Isles represent an ideal field laboratory for tsunami geoscience research. This is due to the widespread preservation of Holocene tsunami sediments in coastal peat deposits. This study uses published accounts of the Holocene Storegga Slide tsunami to illustrate how two different approaches – mapping of tsunami sediments and numerical modelling – produce radically different run‐up heights. The Storegga Slide is one of the world largest submarine slides and took place ca 8150 cal yr bp on the continental slope west of Norway. The tsunami generated by the landslide deposited locally extensive sheets of marine sand and gravel, as well as redeposited clasts of peat across the contemporary land surface. These sediment accumulations have subsequently been buried by peat growth during the Holocene while exposures of the deposits are locally visible in coastal cliff sections. In several areas, the tsunami sediments can be traced upslope and inland within the peat as tapering sediment wedges up to maximum altitudes of between ca 8·1 m and 11·8 m above present sea level. Since reconstructions of palaeo‐sea level for Shetland for ca 8150 cal yr bp suggest an altitude of 20 m below high tide on the day that the tsunami struck, it has been inferred that the minimum tsunami run‐up was locally between 28·1 m (8·1 + 20 m) and 31·8 m (11·8 + 20 m). However, numerical models of the tsunami for Shetland suggest that the wave height may only have reached a highest altitude in the order of +13 m above sea level on the day the tsunami took place. In this paper a description is given of the sedimentary evidence for tsunami run‐up in the Shetland Isles. This is followed by an evaluation of where the palaeo‐sea level was located when the tsunami occurred. Significant differences are highlighted in tsunami inundation estimates between those based on the observed (geological) data and the theoretically‐modelled calculations. This example from the Shetland Isles may have global significance since it exemplifies how two different approaches to the reconstruction of tsunami inundation at the coast can produce radically different results with modelled wave height at the coast being considerably less than the geological estimates of tsunami run‐up.


INTRODUCTION
One of the most problematic areas of tsunami science is how to assess the tsunami risk posed by submarine landslide activity. The issue is well illustrated by the 1998 tsunami in Papua and New Guinea where earthquake modelling could not generate the exceptionally high tsunami run-up observed in the field. Later analysis revealed the existence of a large underwater slump capable of producing such large waves (Syonolakis et al., 2002). Other underwater landslides are known to have occurred in conjunction with offshore earthquakes like the wellknown Grand Banks slide that triggered a tsunami that struck the Burin Peninsula of Newfoundland in November 1929 (Ruffman, 1996) and possibly also a landslide contributed to the Tohoku tsunami in Japan in 2011 (Tappin et al., 2014). The Grand Banks slide was of particular significance since later measurement of the timing of telegraph cable breaks on the seabed led to a reconstruction of the slide and turbidity current velocities across the sea floor (Heezen & Ewing, 1952).
Around 8150 cal yr BP, one of the world's largest submarine slides took place on the continental slope west of Norway (Fig. 1). The slide took place within an area of thick Pleistocene glacigenic deposits, the slide itself is generally thought to have been generated by a strong earthquake in an area located ca 150 km downslope from the Ormen Lange gas field that then developed as a retrogressive slide . The slide produced a huge tsunami that travelled across the Norwegian Sea, Greenland Sea and North Atlantic Ocean. Storegga tsunami deposits have been identified and radiocarbon dated in Scotland, Norway, Faeroe Islands, Shetlands, Denmark and eastern Greenland (Dawson et al., 1988;Bondevik et al., 1997Bondevik et al., , 2003Bondevik et al., , 2005aGrauert et al., 2001;Smith et al., 2004;Wagner et al., 2007;Romundset & Bondevik, 2011;Fruergaard et al., 2015;Rasmussen et al., 2018). The best age estimate is currently 8120 to 8175 cal yr BP that is cited herein as ca 8150 cal yr BP (Bondevik et al., 2012). Separately, radiometric dating of the slide itself yielded similar ages (Haflidason et al., 2005) thus establishing the probability that the tsunami deposits were produced by the tsunami generated by the slide.
The onshore coastal sediments deposited by the Storegga Slide tsunami thus provide a rare opportunity to estimate spatial variations in tsunami run-up and to compare these with wave height values predicted by modelling of the slide and the resulting tsunami. If it is possible to explain the patterns of tsunami by matching the geological field observations with slide modelling, there is an opportunity to apply the results to other continental margins around the world and predict what the patterns of tsunami run-up should be if a future slide were to occur. This particular tsunami has additional interest since the waves passed through an island archipelago and was thus likely to have produced complex patterns of wave inundation. This paper explores these relationships using published accounts of the Storegga tsunami as it affected the Shetland Isles (Fig. 2). This area represents an ideal field laboratory for tsunami geoscience research. This is because of the widespread preservation of Holocene tsunami sediments in coastal peat deposits. This study considers the two different approaches to the reconstruction of tsunami inundation at the coast; numerical modelling and mapping of tsunami deposits.
Where was palaeo-sea level when the Storegga tsunami took place?
The establishment of the position of relative sea level in the Shetland Isles at the time the tsunami took place is a fundamental requirement needed to reconstruct the extent of flood run-up. Unfortunately, however, there is very little sealevel data for Shetland. There are no Lateglacial or Holocene raised beaches on Shetland and the general consensus is that the archipelago was characterized by a history of submergence throughout the Holocene (Hoppe, 1965). The latter author reported from a dredging operation in Lerwick Harbour the presence of tree stumps at a depth of À6Á4 m below high-water mark (HWM) and also in the Sullom Voe area at ca À5Á8 m below HWM. In addition, submerged peat was recovered at Symbister, Whalsay, at between À8Á6 m and À8Á9 m below HWM. A series of five radiocarbon dates on peat and wood obtained by Hoppe (1965) from the above horizon range between 5945 and 6970 14 C yr BP and these were used by Hoppe (1965) to argue that sea level in this part of the Shetland Isles was: "at least about 9 m lower about 6 ka and still at least 9 m lower about 5.5 ka". The radiocarbon dates of Hoppe broadly convert to a range of values between 5990 cal yr BP and 7900 cal yr BP (Bondevik et al., 2005a).
Additional information on former positions of relative sea level during the early and mid-Holocene in Shetland was described by Bondevik et al. (2005a). The latter authors sampled marine sediments from within the Ronas Voe fjord that is partly separated from the North Atlantic by a submerged rock sill between À14 m and À20 m water depth. A radiocarbon date of 7159 to 7520 cal yr BP for the lowermost marine shell within a 3 to 4 m core sequence from these marine sediments provides a minimum (youngest) age for the overtopping of the (now submerged) rock sill by a rising Holocene sea level. In other words, sea level must have been somewhere above À14 to À20 m sometime between ca 7150 cal yr BP and 7520 cal yr BP.
Apart from the above information, other estimates of former positions of relative sea level have been based on: (i) interpretation of Holocene shoreline isobase maps for Scotland; and (ii) geophysical modelling (cf. Bradley et al., 2011). The most relevant shoreline isobase maps for Scotland are by Smith et al. (2012Smith et al. ( , 2017. For example, Smith et al. (2004, fig. 12) produced three isobase maps of land uplift for Scotland since ca 8000 cal yr BP (ca 7000 radiocarbon years BP) using different methods of analysis; although all three of these maps exclude the Shetland Isles, it is clear for two of the maps that the northern Isles (Orkney and Shetland) have experienced at least ca 10 m of subsidence over the last ca 8000 yr. The third map is more equivocal but points to the observation that Shetland did not experience any land uplift over this time interval. Furthermore, a shoreline isobase map for the shoreline produced during the Younger Dryas (the Main Lateglacial Shoreline) points to a relative sea level in the Shetland Isles between ca 11 700 cal yr BP and 12 900 cal yr BP that must have been lower than À35 m (Smith et al., 2004, fig. 11). The estimates of former sea level based on geophysical modelling display a range of sealevel altitudes for the early Holocene in Shetland (Fig. 3). One of the earliest attempts by Lambeck (1995, fig. 6) points to a relative sealevel position for ca 8000 yr BP of between ca À20 to À30 m; modifying an earlier view (Lambeck, 1993) that relative sea level in Shetland at this time most probably was somewhere in excess of 10 to 15 m lower than present. Similar estimates were later derived by Peltier et al. (2002). These estimates of former sea level derived from geophysical modelling and the submerged peat dates were later summarized by Bondevik et al. (2005a) who concluded that sea level at the time of the Storegga Slide was at least 10 to 15 m below the present. The different estimates of the position of relative sea level in the Shetland Isles ca 8000 yr BP is highlighted here since the range of values that have been discussed in the literature each help to define tsunami run-up. In the most extreme scenario where sea level is placed at À30 m at ca 8150 cal yr BP a tsunami deposit reaching up to 5 m above sea level would have to have been associated with a run-up of ca +35 m. Following here as a realistic estimate for the altitude of relative sea level in the Shetland Isles at the time that the Storegga tsunami took place.

Run-up estimates calculated from Storegga tsunami deposits
Sand layers considered to have been deposited by the Storegga tsunami are widespread across Shetland. To date, the majority of sediment exposures have been identified across the northern part of the island chain, in particular North Mainland, Yell and Unst (Fig. 2). In some areas the tsunami deposits are recognizable as sheets of sediment enclosed within peat. Where these occur, for example at Sullom Voe (Bondevik et al., 2003(Bondevik et al., , 2005aSmith et al., 2004) and Whale Firth, Yell (Dawson et al., 2019), the tsunami sediment sheets rise in altitude and become progressively thinner with increasing distance inland. In most respects the sand sheets resemble palaeo-tsunami deposits observed elsewhere in the world (Dawson & Shi, 2000). Storegga tsunami deposits have also been recovered from sediment cores sampled from coastal lakes in Shetland located between 0Á5 to 3Á0 m above present high tide level (Bondevik et al., 2005a). Radiometric dating of these sediments has been extensive and is summarized in Bondevik et al. (2005a) and by Smith et al. (2004). The present study focuses principally on the maximum altitudes of the Shetland tsunami deposits above sea level and what this means in terms of former tsunami run-up, since this factor is critically important to the hydrodynamic modelling of the landslide-generated tsunami.
Storegga tsunami sediment exposures are widespread in coastal peat areas between the Sullom Voe oil and gas terminal and Scatsa airport (Fig. 2). In this area, numerous sediment exposures are visible in peat cuttings while extensive hand-augering has shown that a distinctive tsunami sediment unit is present nearly everywhere; rising in altitude from just below high tide sea level in coastal peat exposures to an altitude of 11Á8 m above sea level where the sheets of sediment have tapered and thinned to such an extent that only individual grains of sediment are visible within the peat (cf. Smith et al., 2004).
A separate suite of Storegga tsunami deposits occurs on the western side of Sullom Voe adjacent to the Houb and Maggie Kettles Loch (Fig. 4). At the Houb the tsunami sediment unit is almost a continuous feature along the shoreline for ca 150 m in coastal peat exposures. However, at Maggie Kettles Loch the tsunami sediment is a visible as a complex deposit of coarse sands and gravels with extensive peat intraclasts to seaward that rises in altitude and becomes progressively thinner inland where it can be traced to 9Á2 m above sea level; it was also traced in hand cores to continue 2Á4 m below high tide in peat underneath the present day beach gravel (Bondevik et al., 2005a).
A third excellently preserved Storegga tsunami sediment occurs on the Atlantic coastline of the island of Yell at Whale Firth (Dawson et al., 2019). Here, tsunami deposits are exposed in coastal outcrops and can be traced upslope as a continuous sediment sheet within peat up to a maximum altitude of 8Á1 m above present sea level. Whilst the existence of such coastal sediment sheets provides convincing evidence for former tsunami deposition, they do not provide evidence of maximum tsunami run-up, since the tsunami almost certainly flooded inland to higher elevations yet did not deposit sediment. This is the case in modern tsunami deposits, for example in Japan 2011 where sand deposits (0Á5 cm thick) were found to around 3 km inland whereas the wave inundated to almost 5 km inland (Abe et al., 2012). In terms of reconstructing tsunami run-up this means that the upper limits of the tsunami sediment sheets simply provide estimates of minimum tsunami run-up.

Storegga tsunami and slide modelling
One of the most difficult aspects of modelling tsunamis generated by submarine slides is converting the model output to run-up values at the coast. The various attempts to model the Storegga tsunami demonstrate this very clearly. In addition to the model of the tsunami generating mechanism, two other models are needed: the first is to model the tsunami as it propagates outward into the open ocean from its source and the second is to model what happens to the tsunami when it approaches and inundates specific coastal areas. These models can be analytical or numerical in nature. In the latter circumstance of a numerical model that attempts to simulate the run-up of the wave onshore, nearshore bathymetric and topographic data are used to create numerical grids of sufficient resolution to provide realistic estimates of how the tsunami waves are deformed within shallow water.
Several attempts have been made to simulate the Storegga tsunami (Harbitz, 1992;Henry & Murty, 1993;Pedersen et al., 1995;Bondevik et al., 2005b;Hill et al., 2014;Lovholt et al., 2017). The first study by Harbitz (1992) made use of a numerical model that had previously been used to simulate waves generated by landslides into fjords in western Norway (cf. Harbitz et al., 1991). The model was then modified to simulate an underwater landslide. It was based on the use of hydrostatic shallow water wave equations to reconstruct wave propagation in the open ocean together with a prescribed slide model that could describe the dynamics of the slide itself. This model of the slide prescribed the shape, run-out distance, maximum velocity, and hence acceleration, all of which affect the wave generated. The second part of the model was based on a comparison of the results from the idealized numerical run-up model using no flux boundary conditions with a separate model for calculating run-up heights on a gentle beach slope (Harbitz, 1992). These calculations were made by Harbitz (1992) for three locations on the Scottish mainland together with one site in eastern Iceland and another in eastern Greenland but no estimates were made for Shetland. It should be noted also that modelling of the local amplification of the tsunami as it approached the nearshore zone was assumed to take place along a linear seabed slope from a chosen point offshore (usually between 70 to 250 m water depth) where the calculated wave height and water depth were known, together with the estimated wave period. In reality, the nearshore seabed almost always is characterized by complex and irregular topography. The run-up calculations for the north of the Scottish mainland suggested that the observed minimum run-up values based on the identification of tsunami deposits was compatible with a Storegga landslide with a velocity of 35 ms À1 , and also consistent with the likely velocity of the 1929 Grand Banks slide and turbidity current (Heezen & Ewing, 1952;Harbitz, 1992). It should be noted also that the modelled slide velocity to 20 ms À1 produced theoretical run-up values that were too low, while a modelled slide velocity of 50 ms À1 produced run-up values that were too high by comparison with the maximum altitudes of observed tsunami deposits in the coastal zone. Henry & Murty (1993) separately attempted to estimate Storegga tsunami run-up by using an estimate of the volume of the Storegga Slide and a simple set of analytical formulae to estimate tsunami amplitude at source; these results indicated a 12 m tsunami amplitude at source. However, although providing run-up values for northern Scotland that are consistent with the geological data for run-up, they did not provide any run-up information for Shetland. The study by Henry & Murty (1993) differs significantly from that of Harbitz (1992) in that by simply using slide volume, they did not simulate the slide itself as part of the calculations of tsunami propagation. It should also be remembered that in both of the above tsunami models use is made of the position of sea level as it exists today and not the very different position it is likely to have been ca 8.15 ka.
More sophisticated numerical models of the Storegga tsunami were later constructed by Bondevik et al. (2005b) and Hill et al. (2014). Bondevik et al. (2005b) used better data on the volume and shape of the Storegga Slide. These authors also looked at how the retrogressive slide motion affected the tsunami generation by performing a simple model where the slide was divided into a series of 167 failing blocks (240 m thick and 600 m long) released one at a time to simulate the retrogressive slide motion (cf. Haugen et al., 2005;Lovholt et al., 2005;Fig. 5). In these models a higher resolution ocean grid was used but the same problem persisted in respect of how best to model the tsunami as it approaches the coastline. Bondevik et al. (2005b) used for such areas a 500 9 500 m grid that was interpolated for certain local areas of western Norway to a 250 9 250 m grid. For Shetland, a mathematical adjustment was employed using Green's law to convert offshore wave height to run-up elevations. The highest estimates of tsunami run-up were made by Bondevik et al. (2005b) who calculated the modelled surface elevations of the tsunami to the north of Shetland as between 7Á1 m and 8Á0 m. Numerical modelling of the tsunami offshore ca 3 km north of Maggie Kettles Loch in a present water depth of 48 m indicates a run-up of 19 to 21 m (Bondevik et al., 2005b). Hill et al. (2014) drew attention to the difficulty of creating an accurate representation of tsunami run-up for areas of complex coastal configuration and bathymetry. These authors attempted to address this problem through the use of multiscale simulation, where grid resolution exhibits spatial variations sufficient to replicate coastal zone topographies but the model did not include inundation of the coastline explicitly. Simulations of the Storegga tsunami by Hill et al. (2014) not only produced a model using present sea level but also made a separate model where an estimate of the palaeo-sea level at ca 8000 cal yr BP was included (Bradley et al., 2011). Similar to Harbitz (1992), the model by Hill et al. (2014) made use of a slide that moved as a prescribed single rigid block; this model using present sea level predicts for Shetland a maximum wave height offshore (referred to as free surface water levels in the model) of ca 8 m. However, when the model is re-run using the palaeo-bathymetry, the free surface level increases by ca 5 m, particularly along parts of the coastline of eastern Shetland. The modelling of Hill et al. (2014) is sophisticated both in terms of the use of a multi-scaling approach and of palaeo-bathymetric data. Its drawback is in the use of an initial block slide rather than the more complex retrogressive slide dynamics envisaged by Bondevik et al. (2005b) and the lack of inundation of the coastline by the waves.
The most recent models are those of Lovholt et al. (2017) who simulated the Storegga Slide both as a retrogressive slide and as a debris flow, and where the run-out of the modelled slide matches the run-out of debris in the slide area. This analysis makes use of the palaeo-bathymetric data of Hill et al. (2014) and is the most sophisticated model to date. The debris-flow slide returned a run-up of 13 m off the coast of Shetland. Still, the models fail to replicate the high tsunami run-up values for Shetland derived from the geological data but they use wave height values offshore rather than simulate the run-up.

DISCUSSION
The research undertaken so far on the Storegga Slide tsunami in the Shetland Isles highlights major discrepancies in the estimates of tsunami run-up at the coast. The geological evidence of tsunami-deposited sand sheets that occur at several locations along the Shetland coastline show that the highest limit of these range between +8Á1 m and +11Á8 m above present sea level. Since it is generally agreed that the upper limit of tsunami sediment deposition is always less than the upper limit of inundation associated with a tsunami, the present authors can be confident that at all of the locations on Shetland where tsunami-deposited sand sheets exist, the maximum run-up was most likely higher than the highest tsunami sediments. Therefore, at Sullom Voe in Shetland where the highest tsunami sediments occur at 11Á8 m above present sea level the inference has to be that tsunami run-up in this area was a minimum of (11Á8 m + 20Á0 m) 31Á8 m. If sea level at ca 8.15 ka was any lower than À20 m, the estimated tsunami run-up would be even higher.
The often marked differences in tsunami runup estimates between those based on geological data and those derived from numerical modelling may be due to a number of factors. The first arises from how the submarine slide is modellednamely whether one uses the single block slide of Harbitz (1992) or Hill et al. (2014) or the retrogressive process of slide failure adopted by Bondevik et al. (2005b). More recently, Lovholt et al. (2017) have addressed this issue in respect of the Storegga Slide and the neighbouring Traenadjupet Slide. These authors have argued that sometimes the use of a retrogressive block model results in lower slide velocities when compared with a single block failure, the latter ultimately resulting in lower tsunami run-up values. More sophisticated slide models may be required to replicate the details of the tsunami which may focus waves on particular areas. Models like those used by Smith et al. (2016) or Abadie et al. (2012) that reproduce deformable slides may provide more insight into the properties of a submarine slide-generated tsunami, but are computationally challenging to perform. However, a more likely cause for the differences in estimated run-up may be due to how tsunami wave behaviour is modelled in nearshore areas. Even the refined multiscale grids used by Hill et al. (2014) may not yet be of sufficient accuracy to capture the ways in which tsunami waves are amplified as they progress into shallow water. Furthermore, none of the models so far has been able to reconstruct the effects of tsunami wave interference on wave amplification and the effects that different tsunami waves in a wave train may have on one another. This requires a fuller understanding of the sea level at the time the tsunami took place as well as any later changes in coastal configuration. For the Shetland archipelago this may be a crucially important issue and one might expect, as elsewhere (for example, Hilo, Hawaii, 1946) that the highest run-up occurs in the lee of individual islands (Mader & Curtis, 1994).
The ongoing research on the nature of Storegga tsunami run-up in the Shetland Isles is hindered by a series of difficulties. First, since the upper limit of tsunami sedimentation at the coast is only a minimum estimate of the real run-up, it is impossible to give a precise estimate of how high the tsunami reached in specific areas. Furthermore, the depositional mechanisms for tsunami sedimentation and how these relate to tsunami wave dynamics are not yet fully understood. Second, it is proving difficult to estimate the position of sea level during the early Holocene at the time the tsunami took place and only a rough estimate can be made. Moreover, in order to reconcile modelling and field data the palaeo-bathymetry/topography needs to account for coastal changes since the Storegga tsunami. This comprises both Holocene relative sea-level changes as well as coastal geomorphological changes, the latter having not been incorporated into any numerical model of the Storegga tsunami. Third, while different numerical models have been very successful in simulating the tsunami in the open ocean, it is proving exceptionally difficult to simulate the amplification of the tsunami as it enters shallow water and ultimately reaches the coast and deposits sedimentthis needs to be a priority for future research. Field investigations of the Storegga tsunami deposits in Shetland together with the results of various modelling experiments have provided much valuable information to date on how and why the tsunami took place. However, for the various difficulties in interpretation there are at present no satisfactory answers.

CONCLUSIONS
There are very clear differences in the Shetland Isles between the highest elevations of the Storegga Slide tsunami run-up deposits and the highest run-up estimates generated from every numerical simulation that has so far been developed. In the three areas where the field evidence for tsunami sediment sheets is very clear (Sullom Voe, Whale Firth and Maggie Kettles Loch) the highest altitudes for the deposits range between 8Á1 m and 11Á8 m above present sea level. Bearing in mind that the highest tsunami deposits onshore represent minimum values of run-up and also that this study has used a realistic estimate (À20 m) that regional sea level in the Shetland Isles could possibly have been ca 8150 cal yr BP, minimum tsunami run-up for these three locations ranging between 28Á1 m and 31Á8 m has been derived here. It is almost certain that the tsunami run-up values are higher than these values but by how much more it is presently impossible to determine. It should be noted that the position of relative sea level ca 8 ka cited here as À20 m may have been as low as ca À30 m (Fig. 3A).
Most of the numerical models that have been developed for the tsunami fail to predict such high run-up values. It is suspected here that this problem occurs because the models are not able to replicate the local amplification of the tsunami waves as they enter shallow water at the coast. Furthermore, the models that exist at present are not sufficiently powerful to model the diffraction, refraction, reflection and interference of different waves within the tsunami wave train as it approached the palaeo-coastline. This is challenging computationally due to the need to model large regional areas at the same time as focussing on smaller-scale regions. Coupling studies such as that done by Bondevik et al. (2005b) may prove fruitful in gaining a fuller understanding of wave dynamics in the Shetlands. This may prove to be crucial in respect of a large tsunami reaching an island archipelago where current knowledge of tsunami hydrodynamics highlights the fact that the areas of highest run-up may be in the lee of individual islands. The example from the Shetland Isles described here may have global significance because it exemplifies how two different approaches to the reconstruction of tsunami inundation at the coast can produce radically different results, with modelled wave height at the coast being considerably less than the geological estimates of tsunami run-up.