Sedimentary evidence of the Late Holocene tsunami in the Shetland Islands (UK) at Loch Flugarth, northern Mainland

Tsunami deposits around the North Sea basin are needed to assess the long‐term hazard of tsunamis. Here, we present sedimentary evidence of the youngest tsunami on the Shetland Islands from Loch Flugarth, a coastal lake on northern Mainland. Three gravity cores show organic‐rich background sedimentation with many sub‐centimetre‐scale sand layers, reflecting recurring storm overwash and a sediment source limited to the active beach and uppermost subtidal zone. A basal 13‐cm‐thick sand layer, dated to 426–787 cal. a CE based on 14C, 137Cs and Bayesian age–depth modelling, was found in all cores. High‐resolution grain‐size analysis identified four normally graded or massive sublayers with inversely graded traction carpets at the base of two sublayers. A thin organic‐rich ‘mud’ drape and a ‘mud’ cap cover the two uppermost sublayers, which also contain small rip‐up clasts. Grain‐size distributions show a difference between the basal sand layer and the coarser and better sorted storm layers above. Multivariate statistical analysis of X‐ray fluorescence core scanning data also distinguishes both sand units: Zr, Fe and Ti dominate the thick basal sand, while the thin storm layers are high in K and Si. Enriched Zr and Ti in the basal sand layer, in combination with increased magnetic susceptibility, may be related to higher heavy mineral content reflecting an additional marine sediment source below the storm‐wave base that is activated by a tsunami. Based on reinterpretation of chronological data from two different published sites and the chronostratigraphy of the present study, the tsunami seems to date to c. 1400 cal. a BP. Although the source of the tsunami remains unclear, the lack of evidence for this event outside of the Shetland Islands suggests that it had a local source and was smaller than the older Storegga tsunami (8.15 cal. ka BP), which affected most of the North Sea basin.

Tsunami deposits are an integral component used in the assessment of the long-term hazard of tsunami events, in particular in coastal regions with a short and fragmented historical record (Dawson & Shi 2000;Switzer & Jones 2008;Gonz alez et al. 2009;Weiss & Bourgeois 2012;Engel et al. 2016).Although a range of indicative sedimentary, geochemical and micropalaeontological signatures have been established (Dawson & Shi 2000;Switzer & Jones 2008;Engel et al. 2016; Chagu e 2020; Spiske 2020), the identification of tsunami deposits in the coastal sedimentary record poses challenges at the site scale.Tsunami deposits may vary considerably in structure and thickness depending on the hydrodynamics of the tsunami, the availability and nature of source sediments, the coastal bathymetry and topography, the depositional setting (i.e.onshore, lacustrine, offshore) and potential post-depositional changes.
Despite an apparent scarcity of tsunamis in the shallow North Sea basin, the Shetland Islands have become an important field laboratory in which to study tsunami deposits since the early 1990s (Smith 1993;Bondevik et al. 2003Bondevik et al. , 2005;;Dawson et al. 2006Dawson et al. , 2020b;;Costa et al. 2015;Cascalho et al. 2016;Buck & Bristow 2020), with summaries of the findings presented in Smith et al. (2004Smith et al. ( , 2019)), Long (2015), Costa et al. (2021) and Bondevik (2022).The original studies on Shetland have played a substantial role in the advancement of tsunami sedimentology and wider tsunami geoscience.Most studies since have focused on the abundant sedimentary evidence of the tsunami triggered by the Storegga submarine slides at the Norwegian shelf margin dated to 8175-8120 cal. a BP (Bondevik et al. 2012).The Storegga deposits are usually preserved in thick coastal sections of blanket peat along many of the Shetland fjords, locally called 'voes'.They mostly consist of normally graded sand up to 40 cm thick, with large rip-up clasts, gravel components and pieces of wood embedded in the sand (Bondevik et al. 2003), as well as increased concentration of heavy minerals at the base (Costa et al. 2015;Cascalho et al. 2016).The deposits thin and fine inland, with sharp, in some cases erosional lower boundaries.The maximum run-up of the Storegga tsunami on the Shetland Islands is estimated to be up to 20-25 m above contemporary sea level (Bondevik et al. 2003(Bondevik et al. , 2005;;Dawson et al. 2020b).Deposits of the Storegga tsunami are also found in coastal lakes, consisting of normally graded coarse sand and fine gravel as well as lake-mud clasts and organic macro-remains above an erosive base.There is further evidence of a c. 5500 cal. a BP tsunami deposit found at Garth Loch at the east coast of Mainland (Bondevik et al. 2005), which was possibly also identified in western Norway (Romundset et al. 2015).The youngest tsunami deposit identified so far on the Shetland Islands was dated to 1500 cal. a BP and occurs as a thin inland-fining sand layer at Dury Voe (Bondevik et al. 2005) and Basta Voe (Dawson et al. 2006).Despite the solid evidence, Long (2015) considers the event 'uncertain' and points to possible alternative depositional processes (e.g.aeolian).The potential trigger mechanism for the late Holocene tsunami, however, is unknown (Bondevik et al. 2005;Dawson et al. 2006;Ballantyne et al. 2018).
Here, we present the sedimentary record of a Late Holocene high-energy flooding event from coastal Loch Flugarth, on the Shetland Islands.We discuss its origin based on a range of sedimentary criteria typically observed in tsunami deposits and the deposit's possible relation to the other inferred tsunamis on Shetland.Our objective is to add to the palaeotsunami record on Shetland and to improve the chronology.Finally, we aim at stimulating further investigations of tsunami trigger mechanisms to improve coastal hazard assessment in the North Sea basin.

Physical setting
The Shetland Islands (1.3°W, 60.4°N) comprise a group of more than 100 islands, located at the northern margin of the North Sea, about 165 km NE of the Scottish mainland (Fig. 1A) (Mykura et al. 1976;Bennett et al. 1992).The geology of the Shetland Islands is dominated by metamorphic and metasedimentary rocks of the Caledonian Orogeny.The archipelago is tectonically divided by the major north-south running fault system of the Walls Boundary Fault (Mykura et al. 1976;Gillen 2003), which represents a continuation of the Great Glen Fault (GGF) (Fig. 1A).The GGF extends in southwestern direction all the way through the Scottish mainland (Pringle 1970;Gillen 2003).

Loch Flugarth
The study site of Loch Flugarth is located in the North Roe area (northern part of Mainland; Fig. 1B) and is embedded into the Sand Voe Group, comprising thick banded hornblendic and granulithic gneisses with thick lenses of amphibolite and metagabbro near the top of the Sand Voe Group (Mykura et al. 1976).The lake is crossed by the Flugarth Fault running from SW to NE (Pringle 1970).
The freshwater lake of Loch Flugarth covers an area of 0.16 km 2 .It is very shallow with a maximum water depth of ~2.4 m (Murray & Pullar 1908) and no thermal stratification (Murray & Pullar 1908;Bennett et al. 1992).Such polymictic lakes with a non-varved record, a very simple and shallow bathymetry and limited fetch are ideal to preserve allochthonous event layers as disturbing internal processes such as hyperpycnal flows, interflows or erosion from the flanks are insignificant (Schillereff et al. 2014).
The lake level is ~2 m above mean sea level inside the very shallow and usually calm Sand Voe bay, from which Loch Flugarth is separated by a ~70 m wide sandy barrier (Chapelhow 1965).The barrier, which is crossed by an artificial cut, forms a flat beach with mostly sands and some gravel (Flinn 1974); its elevation is ~2 m above high tide.Stabilised aeolian dunes and till framing the beach show signs of recent erosion, whereas only the largest storm surges manage to overtop the barrier.The barrier seems to have formed inside Sand Voe bay first as a spit fed by the longshore drift, eventually forming a closed barrier.The area behind the barrier was then cut off and turned into a freshwater lake gradually filling with sediment (Flinn 1964(Flinn , 1974) ) from organic detritus, lake-internal bioproduction, potentially aeolian input and recurring coastal overwash processes (Hess et al. 2023b).The only creek feeding Loch Flugarth is the Beorgs of Skelberry in the south, the organic-rich sediment input of which creates a flat, delta-type bathymetry (Murray & Pullar 1908).

Sea-level changes
Changing relative sea levels (r.s.l.) control the onshore accommodation space and sensitivity of coastal lakes to receive sediments during overwash by tsunamis or storm surges and waves (Liu 2004;Szczuci nski 2020).Thus far, Holocene sea-level index points are scarce in Shetland.Based on very few dated basal peats and a limiting shell from Ronas Voe to the SW of the study site, the r.s.l. is assumed to have risen at a high rate until c. 7000 to 6000 cal. a BP and subsequently slowed down (Hoppe et al. 1965;Bondevik et al. 2005;Smith et al. 2019).For the last 1500 years, r.s.l.rise is assumed to have been minimal (Bondevik et al. 2005;Dawson et al. 2006) with tide gauge measurements from Lerwick showing no trend in r.s.l.since c. 1950 CE (Wahl et al. 2013).The tidal range is between 1 m (neap) and ~2 m (spring) (data from nearby Sullom Voe, Halliday 2011).Therefore, we infer a relatively stable geomorphic environment for overwash across the coastal barrier into Loch Flugarth during extreme-wave events of the last c.1500 years.Under such conditions, small coastal lakes are excellent traps for tsunami deposits with good preservation potential (Kelsey et al. 2005;Kempf et al. 2017;Dawson et al. 2020a).

Field methods
In the central northern part of Loch Flugarth at a water depth of ~2.0 m, three sediment cores (FLUG 2-4) with lengths of up to 91.7 cm were taken at a distance of <30 m to each other using a zodiac and a UWITEC gravity corer (Figs 1D, S1, S2).They were kept in transparent PVC liners (60 mm Ø; Fig. S2).Core FLUG 1 was taken using a Russian chamber corer that was opened at a depth of 200 cm below the lake bottom to generate an undisturbed stratigraphic sample at a core depth of 200-250 cm (Fig. S4).Additionally, two surface samples of modern sedimentary environments were taken to identify sediment sources.In this study, we focus on the basal part of cores FLUG 2 and 3.The upper parts of these cores were analysed for storm overwash by Hess et al. (2023b).

Sedimentary analyses
After splitting and describing the cores at the laboratory of the Geological Survey of Belgium at Brussels, a Geotek multi-sensor core logger was used at the Renard Centre of Marine Geology, Ghent University, Belgium.High-resolution photographs were taken using the Geoscan IV line-scan camera.Bulk density (c-ray attenuation densometer) and magnetic susceptibility (Bartington point sensor MS2E) were measured in steps of 2 mm.X-ray computer tomography (CT) scans were generated at Ghent University Hospital using a Siemens Somatom Definition Flash Medical X-ray CT scanner.The voxel sizes of the image stacks are 0.15 9 0.15 9 0.3 mm.The stacked core images were analysed qualitatively using the Fiji software package (Schindelin et al. 2012).
For grain-size analysis of the bottom sand unit in FLUG 3, a Fritsch Analysette 22 NeXT Nano laser diffractometer with a measurable range of 0.1-3800 lm was used at the Laboratory for Geomorphology and Geoecology, Heidelberg University, Germany.The thickness of each sample was 2 mmone-fifth of what is recommended for tsunami deposits (Spiske 2020).Samples were dried at 105 °C and carefully pestled.Organic matter was dissolved by adding H 2 O 2 (30%).Aggregates were dispersed using 0.1 mol Na 4 P 2 O 7 (44.6 g L À1 ).GRADISTAT v9.1 (Blott & Pye 2001) was used to calculate univariate statistical measures following Folk & Ward (1957).
Elemental concentrations were measured with an Avaatech (GEN-4) X-ray fluorescence core scanner at the Institute of Earth Sciences, Heidelberg University.The core scanner uses an OXFORD 'Neptune 5200' series 100 W X-ray source with a Rh anode.The slit size was 5.0 mm down-core and 10.0 mm cross-core.Two different energy levels were applied for measurement: at 10 kV with a current of 500 lA without a filter and a dwell time of 10 s, and at 30 kV beam with a 1500 lA current, a Pd-thick filter and dwell time of 10 s.The data were processed using the bAxil software (Brightspec 2015).The elements were normalised by dividing the intensity counts of each element by the X-ray fluorescence total during each run excluding coherent and incoherent scattering (Kern et al. 2019).
Principal component analysis (PCA) was first carried out on 13 centred log-ratio transformed elements (Al, Br, Ca, Fe, K, Mn, Mo, P, Rb, S, Si, Sr and Zr) and subsequently on indicative element ratios.The first five principal components (PCs) from the second PCA were used for the hierarchical clustering of principal components (HCPC) to identify the most significant variables and to identify potential differences between types of sand layers in the cores.The PCA and the HCPC were performed in R studio (Kassambara 2017) by applying the R packages FactoMineR (v2.4) for computing and factoextra (v1.0.7) for visualisation.Significance tests were applied to analyse whether the differences between the clusters are statistically significant.A Kruskal-Wallis test (R package vegan v2.5-7) and a Dunn's post-hoc test (R package dunn.testv1.3.5)identified characteristic variables controlling the clusters.The output resulted in Bonferroni-corrected p-values with a > 0.05.
A total of nine samples from selected core depths representing the different facies were processed for palynological analysis following standard techniques, including: freeze-drying, weighing (range of sample dry weights 0.1-0.7 g), spiking with Lycopodium spores (batch no.050220211), treatment with HCl, NaOH and HF, heavy liquid separation, sieving (10 lm) and mounting of the residues on glass slides using glycerin jelly.All samples were screened to identify the most common palynomorphs.

Chronology
A combination of 14 C and 137 Cs dating was used to establish a Bayesian age-depth model for the entire stratigraphy of FLUG 2 and 3. Samples for 14 C dating were measured by accelerator mass spectrometry at the Pozna n Radiocarbon Laboratory, Poland (Goslar et al. 2004).Calibration of the 14 C data was carried out using the IntCal20 data set (Reimer et al. 2020).Details on 137 Cs dating are compiled in Data S1.For Bayesian age-depth modelling rBacon (Blaauw & Christen 2011) was applied.Variable sedimentation rates or event layers were considered in the model according to Blaauw & Christen (2011) (Table S2).Late Holocene 14 C dates of tsunami stratigraphies in the Shetland Islands presented in Bondevik et al. (2005) were recalibrated using Calib software (v8.2) (Stuiver & Reimer 1993) and the IntCal20 data set (Reimer et al. 2020) (Table S3).Data on the bulk density, X-ray fluorescence and CT scans and the agedepth model were first published in Hess et al. (2023b).
The modern surface samples representing potential source environments for ex-situ deposits inside the lake were taken from the foot of the barrier at the lake margin (FLUG-M1) and from the wet beach (FLUG-M2) (Fig. 1D).Both are very similar, consist of pure, relatively well-sorted sand and can be considered representative of the major part of the barrier complex.
The thick basal sand and the remaining 40 thin sand layers above, which were identified based on the CTscans (Fig. 2), share a similar sand source.However, the multiproxy data reveal differences that go beyond the mere thickness.Compared with the thin sand layers above, the basal sand is finer and less sorted, and free of carbonate, and the pore space is partially filled by a mud matrix.
The magnetic susceptibility is highest in the basal sand, much higher than in the thinner sand layers.Elemental ratios of Sr/Br (Fig. S3) and Zr/K reach maximum values in the basal sand section and are up to almost an orderof magnitude higher than in the thin sand layers above.In contrast, S/Ti values, which peak in facies 1 (muddy peat), are lowest in the basal sand, much lower than in the thin sand layers L1-L40 (Fig. 2).
These differences are corroborated by the PCA and the HCPC.Initially, element counts were used for the PCA.The first two PCs explain 47.2% (PC1) and 14.0% (PC2) of the variance of the whole data set, respectively (Fig. S5).PC1 is driven by the variance of Al, Br, K, Rb, S, Si and Zr, whereas PC2 is mainly driven by Ti.These results were used to establish ratios with conservative elements (e.g.Ti, Fe, Si, K) as the denominator, mostly Ti as it is neither contributed by internal biogenic processes nor influenced by redox conditions or diagenetic overprinting (Chagu e 2020).The PCA on ratios shows four separate clusters and a clear distinction of the basal sand (cluster 1) from the remaining sand layers (cluster 2) and the organic-rich background facies (clusters 3 and 4) (Fig. S6).While the basal sand shows high ratios of Zr/K, Sr/Br and Zr/Rb, the thin sand layers show high ratios of K/Ti, Si/Zr and Si/Ti.The HCPC based on the eigenvalues reveals five dimensions and also four clusters, very similar to those of the PCA.The sands of facies 4 are again clearly separated into clusters 1 and 2 (Fig. 3, Table 1) with Zr/K and Sr/Br showing the strongest significance for the basal sands (cluster 1).In contrast, the upper thin sand layers (cluster 2) are mostly driven by Si/Ti and K/Ti (Fig. 4).
An initial screening of pollen content across core FLUG 3 also revealed differences between the basal sand and the remaining core: the seven samples between 3 and 73 cm b.s. are all rich in terrestrial pollen and spores.They are dominated by Ericaceae (heath), Poaceae (grasses) and Cyperaceae (sedges) pollen representing the heath and grasslands around Loch Flugarth.These samples also contain various pollen grains from herbs and fern spores, as well as sporadic arboreal pollen such as Betula (birch), Alnus (alder), Pinus (pine) and Salix (willow).In contrast, the two samples from the basal sand (81 and 86 cm b.s.) contain only very few reworked pollen grains of poor preservation.

Sublayers of the basal sand
The grain-size analysis allows for an internal differentiation of the basal sand into four sublayers (Fig. 5).From bottom to top, sublayer 1 (91.7-88.5 cm b.s.) has a lower section of unimodal, moderately to poorly sorted massive sand, and highest values of Zr/K and bulk density.The upper part is normally graded, shifting to bimodal (sand and minor silt component) and trimodal (sand and minor silt and clay components) grain-size distributions (Fig. 7), which are very poorly to extremely poorly sorted.The Zr/Kvalues and bulk density decrease in the upper part.
Sublayer 2 (88.5-82.7 cm b.s.) has a lower inversely graded section changing from a mixture of sand and silt to pure sand, with better sorting towards the top and slightly increasing bulk density as well as stable magnetic susceptibility and Zr/K values.The upper part is again normally graded with a distinct change from silty sand to sandy silt at 83.9 cm b.s. and the highest values for magnetic susceptibility of the entire succession.The uppermost part of sublayer 2 is less sorted and shows a decrease in Zr/K and an initial doubling of S/Ti ratios.This might be associated with millimetre-scale muddy peat clasts embedded in the silty to sandy matrix, the largest of which is >1 cm in vertical extent (Figs 5, S8).
There is a sharp contact to sublayer 3 (82.7-80.6 cm b.s.) with a very thin inversely graded base followed by a normally graded top section with increasing concentration of organic matter and small mud clasts.The Zr/K decreases in the entire sublayer, along with bulk density and magnetic susceptibility, while the S/Ti levels remain elevated.
Sublayer 4 (80.6-78.6 cm b.s.) of massive sand has an erosive lower boundary, elevated S/Ti values and millimetre-scale muddy peat clasts.It shows a distinct peak in magnetic susceptibility and bulk density

Thick basal sand
Tsunami deposition associated with a weak increase in Zr/K values, and is better sorted than sublayer 3. Its upper boundary to facies 3 is sharp, where Br/Ti (Fig. S3) and S/Ti rise strongly, while the magnetic susceptibility, bulk density and Zr/K values decrease (Fig. 5).

Age of the basal sand
Age-depth modelling based on 14 C and 137 Cs reveals a chronostratigraphy of c. 1450 years.The four 137 Cs dates (see age model in Data S1) and the surface date were integrated together with the 14 C data into the main agedepth model for the entire core.Two outliers (FLUG 3-30-31 and FLUG 3 14C-4; Fig. 6, Table 2) both creating significant age inversions in the stratigraphy were not considered for further analyses.After removing all depositional events (in total 19 sand layers, L ≥ 2 mm; Table S2) by identifying their upper and lower boundaries in the age-depth model and subsequent reintegration at their real depth (cf.(Fig. 6).Further details on the construction of the age model can be found in Data S1.

Sediment sources
The three types of organic-rich facies mainly represent lake-internal sedimentation with a contribution of organic matter through surface discharge from the peat-bog catchment.The highest Br and S in facies 1 (cluster 4), which are best preserved in fine-grained organic matter (Engel et al. 2012 The thin sand layers of facies 4 (L1-L40) have been interpreted as deposits of individual overwash events originating from the dune, beach and uppermost subtidal area of Sand Voe during major storm events.This interpretation is mostly based on their similarity with local beach sands, sharp upper and lower contacts and occasional gravel components (Hess et al. 2023b).The combination of tides, storm surge and wave run-up exceeding the threshold height of the barrier displace sand from the beach and the barrier and wash it into the coastal lake.They may create washover fans proximal to the barrier, mainly through bedload sedimentation, which extend to thin centimetre-to millimetre-scale sand sheets inside the coastal lake (Liu 2004;Chaumillon et al. 2017).The sediment source area for the thin sand layers is confirmed by the large overlap of endmember 1 of the endmember modelling analysis of grain-size distributions of the entire core by Hess et al. (2023b).Endmember 1 represents the thin layers of sandy storm deposits and fully overlaps with the grain-size distribution of the modern beach sand and the subaerial backside of the barrier (Fig. 7).Endmember modelling decomposes grain-size distributions of a data set by calculating a low number of synthetic distributions best representing the variability of the data set.This helps to distinguish and group similar distributions and to unmix and quantify different sediment sources represented in individual samples (Dietze & Dietze 2019).
Facies 4 is separated into two distinct clusters, i.e. the thin storm sand layers and the thick basal sand layer.The basal sand is dominated by high values of Zr/K, Fe/Ti and Sr/Rb.Zr/K most likely reflecting the concentration of zircon as one of the highest-density minerals (Davies et al. 2015), which, despite generally low percentages, has been shown to be enriched in tsunami deposits (Costa et al. 2015; Chagu e 2020).Zircon grains require highenergy flows to generate sufficient shear stress to become mobilised (Cuven et al. 2013).Furthermore, the very low S/Ti ratio (Fig. 4) may be driven by high Ti concentrations associated with heavy mineral content (cf.Cuven et al. 2013; Chagu e 2020).In contrast, the storm deposits are dominated by K and Si, i.e. mainly quartz and feldspar varieties (Chagu e 2020) as in the beach sand.This is confirmed by the HCPC (Fig. 3) as well as the Kruskal-Wallis test and Dunn's post-hoc test (Fig. 4).

Evidence for tsunamigenic origin of the basal sand
The basal sand stands out against the background sediments in terms of bulk density and magnetic susceptibility, which is caused by high-density ferromagnetic iron oxides or heavy minerals, and is a typical signature of tsunami deposits in lakes (Wagner et al. 2007;Kempf et al. 2017).Suspension grading, PC1 is mainly driven by Br/Ti and S/Ti and, thus, by organic matter.It inversely correlates with heavy mineral content.This may be a function of sediment source from distal (shallow marine to coastal) with low values to increasingly proximal (beach ?catchment ?internal bioproduction).PC2 is mainly driven by Si/Ti and Si/Zr (Fig. S6), reflecting the beach (and dune) sand.
which is observed in sublayers 1-3 (Fig. 5), is a common sequential pattern in tsunami deposits (Bondevik et al. 1997;Dawson & Shi 2000).It forms when a tsunami loaded with suspended sediment impacts a coastal lake or a coastal plain still inundated from the previous wave.Grains settle as a function of their size, shape and density under Stoke's law of hydraulic equivalents, with smaller and lighter particles increasing towards the top (Dawson & Shi 2000;Spiske 2020).Therefore, the heavy mineral content, in the current case reflected by Zr/K, tends to be higher in the lower parts of sublayers that are only partially formed by suspension grading (cf.Bahlburg & Weiss 2007;Moore et al. 2011;Cascalho et al. 2016;Spiske 2020).The increased heavy mineral content in combination with the bi-to trimodal pattern of finer sand grains and mud in samples from the upper part of the normally graded sublayers (Figs 5, 7) indicates a larger source area than the thin sand layers above.This source area probably comprises the dune, the beach, Sand Voe bay and marine areas outside of Sand Voe bay, at water depths out of reach for storm waves.This pattern of the source area is typical of tsunami deposits when compared with storm deposits (Switzer & Jones 2008;Engel et al. 2016;Spiske 2020).
In between sublayers, sharp and undulating erosional boundaries are common (Cuven et al. 2013;Spiske 2020), such as at the base of sublayers 3 and 4 (Fig. 5).Tsunamis erode underlying fine-grained and/or organic-rich substrate onshore during both run-up and backwash.The cohesive clasts become entrained in the turbulent flow and embedded into the sandy tsunami deposit as rip-up clasts (Cuven et al. 2013;Spiske 2020), which have been reported from tsunami deposits in various coastal lake archives (Bondevik et al. 1997;Kelsey et al. 2005;Wagner et al. 2007;Kempf et al. 2015Kempf et al. , 2017)).
Sublayer 2 shows a distinct inversely graded section, which may relate to the process of kinetic sieving where grain-grain collisions in the basal part of the flow cause finer particles to trickle down and settle before the larger particles (Sohn 1997).Such traction carpets indicate high shear stress from the suspension-laden flow and have been identified in several tsunami deposits (Moore et al. 2011;Falvard & Paris 2017).However, an upward reduction in bulk density with increasingly coarser particles as observed by Moore et al. (2011) is not reflected by the present data.
A distinct mud cap on top of the entire tsunami succession, such as observed in an organic-poor, gyttjatype lake (Kempf et al. 2015(Kempf et al. , 2017) ) or clay-rich marsh deposit (Cuven et al. 2013), cannot unequivocally be inferred from the photographs and CT scans (Fig. 2).However, there is an increasing concentration of finer particles and organic matter in the upper part of sublayer 3, possibly reflecting the waning stage between two waves, a pattern that has been found in other tsunami deposits of shallow coastal lakes as well (Bondevik et al. 1997;Bondevik 2022).The low concentration and poor preservation of pollen in the basal sand compared with all other facies of the core are in agreement with observations from tsunami deposits in the region (Smith et al. 2004) and elsewhere (Chagu e-Goff et al. 2012).
The characteristics of the basal sand section overlap with typical criteria of tsunami erosion and deposition in the proximal coastal lake environment sensu Kempf et al. (2017).Alternative mechanisms can be excluded, as we can clearly distinguish them from the site-specific pattern of the overlying storm deposits, while the sorting is too poor for aeolian input.Furthermore, we exclude major changes to the geomorphic framework of the sedimentary archive of Loch Flugarth, as r.s.l., which is the most crucial controlling factor, is assumed to have not changed notably during the last 1500 years (Dawson et al. 2020b).
The investigated gravity cores may not reach the base of the candidate tsunami deposit.However, this affects neither the interpretation presented above nor the chronological estimates discussed in the following section.The finding of brownish organic-rich (comparable to facies 1-3) layers and some greyish, more minerogenic layers at a depth of 250-200 cm below the lake bottom in FLUG 1 (Fig. S4) shows that the record may extend beyond the Late Holocene.Prior work at the site indicates two distinctive units deeper down in the record possibly corroborating with earlier highenergy inundation to the lake basin (Sue Dawson, unpubl.data).Thus, using more advanced coring techniques, a more complete Holocene record and additional event deposits may be accessed and investigated at Loch Flugarth.

Timing and trigger of the tsunami
The age range for the basal sand of 426-787 cal. a CE or 1524-1163 cal. a BP, respectively (536-701 cal. a CE or 1414-1249 cal. a BP with 1r error), is based on the age-depth model and considered to represent the time window in which the inferred tsunami occurred.Age overestimation owing to excessive reworking of the 14 C-dated material used for the age-depth model is assumed to be negligible, as outliers were manually eliminated before running the age-depth model (Data S1, Fig. 6, Table 2) (Hess et al. 2023b).
There are tsunami deposits found at two locations on the Shetland Islands which overlap with this time range.At Dury Voe, east Mainland, a thin sand layer located within thick coastal peat was traced in outcrops from the inner fjord up to some 400 m inland and up to an elevation of 5.6 m above high tide.Two 14 C dates were retrieved 7-8 cm (2303-1927 cal. a BP) and 1-2 cm (1818-1520 cal. a BP) below the sand layer.One 14 C dating was generated 1-2 cm above the sand layer (1507-1287 cal. a BP) (Table S3) (Bondevik et al. 2005).These data are maximum ages as the dated objects could have been subject to some limited reworking.Also, the lower ages pre-date the event by an unknown amount of time as it is unclear how much of the peat stratigraphy was eroded by the tsunami.The location is proximal to the shoreline where shear stress by tsunami flow may still have been high.We argue that the maximum age generated right on top of the sand layer is more indicative for the timing of the Dury Voe tsunami event and may suggest an age younger than the 1500 cal. a BP estimated by Bondevik et al. (2005).
At the inlet of Basta Voe, northeast coast of Yell, three thin sand layers were identified in peat outcrops close to the coast, the uppermost of which continuously extends for ~2 km inland (Dawson et al. 2006).It can also be traced in ground-penetrating radar measurements (Buck Table 2. 14 C data of samples from Loch Flugarth (Hess et al. 2023b).Calibrated ages with an asterisk indicate deposition after 1950 CE.Those in brackets were interpreted as outliers (too old) and not considered for the age-depth model. 14C data were calibrated using the IntCal20 data set (Reimer et al. 2020) & Bristow 2020).The sand layer is bracketed by two overlapping 14 C dates: 0-1 cm below the sand layer, 1546-1345 cal. a BP; 5 cm above the sand layer, 1517-1303 cal. a BP (Table S3) (Bondevik et al. 2005;Dawson et al. 2006).Based on the overlap, the age range probably covers the event.Chronological data from both sites point to an age of the tsunami around 1400 cal. a BP (or even younger) rather than 1500 cal. a BP, considering that reworking of dated material may also have played a role.This chronological estimate is in accordance with the chronological interpretation of the basal sand in FLUG 2 and 3 (Fig. 8), which possibly represents the Dury Voe tsunami sensu Bondevik et al. (2005).
A source for the Dury Voe tsunami has not yet been identified.A local coastal or submarine landslide has been discussed (Bondevik et al. 2005;Dawson et al. 2006;Ballantyne et al. 2018), but without any physical evidence from seafloor mapping (Long 2015).Based on their limited size, none of the post-Storegga slides known north of the Shetlands seems capable of triggering a tsunami (Haflidason et al. 2005;Smith et al. 2019).contain the record of 'An earthquake in October'.Both events are ascribed to the year 720 CE, indicating a major seismic event triggering a tsunami that impacted the Outer Hebrides and potentially other coastal regions of the northern British Isles.Mac Conamhna (2023) associates this event with seismicity along the GGF running across the Scottish mainland from the southern Hebrides up to the Shetland Islands in a SW-NE direction, and draws a correlation with the deposits of the Dury Voe tsunami on Shetland.However, Fig. 8 shows that maximum ages for the event beds at Dury Voe, Basta Voe and Loch Flugarth appear too old.There is a small overlap with the age range of the Flugarth deposit from this study, although the 720 CE archive entry is more than 160 (2r error) or 126 (1r error) years younger than the median age from Flugarth.A correlation can only be established under the assumption that the dated objects at all three sites are systematically older than their sedimentary contexts, which is considered unlikely.

Comparison with the historical record of tsunamis
There are diverging conclusions on the earthquake record of the GGF (Musson 2007;Piccardi 2014), which is mostly driven by glacio-isostatic rebound (Davenport et al. 1989;Ringrose 1989).The two strongest historical earthquakes on Scottish territory with proper documentation, both of which could be associated with the GGF, have estimated magnitudes of 5 (1901 CE) and 5.1 (1816 CE), respectively (Piccardi 2014), with no reported tsunamis.Onshore palaeoseismological data indicate maximum possible magnitudes of 6.0-7.5 following deglaciation (13 000-6000 cal. a BP) and of 6.5 afterwards (Davenport et al. 1989).Based on these data, it is highly speculative whether a submarine segment of the GGF had the ability to generate a major earthquake and tsunami that impacted coastal sections spanning from the Outer Hebrides up to the Shetland Islands and generate a tsunami bed at several sites with both northern and eastern exposure across the Shetland archipelago.Thus, the trigger mechanisms of the Dury Voe tsunami and the seismic and coastal flooding event recorded near Iona still remain elusive based on the few existing records.

Conclusions
Sediment cores from Loch Flugarth on northern Mainland, Shetland Islands, show a clear separation between organic-rich mud of varying composition and thin allochthonous sand layers.The allochthonous sand occurs in thin sub-centimetre layers throughout the core and results from storm overwash.These thin layers differ from the thick basal sand in all cores in terms of grain-size distribution, sediment structure and chemical composition.High values of Zr/K and magnetic susceptibility as well as low S/Ti values point to enriched heavy minerals.Poorer sorting, partially with bi-and trimodality, as well as suspension grading in sublayers, some of them also showing traction carpets and rip-up clasts, indicate deposition by a tsunami.The age-depth model and correlation with tsunami deposits of similar age in the Shetland Islands suggests an age around 1400 cal. a BP (c.536-701 cal. a CE; 1r error).Correlation with a historical, allegedly seismic tsunami observed on the Outer Hebrides in 720 CE seems questionable based on the chronological evidence.
While the present study adds to the growing evidence of the Late Holocene Dury Voe tsunami, it prompts the need for further research on the spatial distribution and impact of tsunamis in the North Sea region, and the possible triggers as a basis for future hazard assessment.From a more local perspective, this study corroborates the assumption that the site represents an archive of event deposits beyond the last 1500 years.Loch Flugarth preserves both tsunami and storm deposits that can be unequivocally distinguished based on sedimentological and geochemical criteria.Therefore, we advocate for exploration of the deeper sediment record of Loch Flugarth to improve our chronological estimate and correlation of the Dury Voe tsunami and possibly other regional tsunamis such as the Garth tsunami (c.5500 cal. a BP) and to also extend the palaeotempestological record of the northern North Sea.Table S1.Results of the mathematical 137 Cs modelling (Hess et al. 2023b).b.s.= below surface.
Table S2.The 19 most pronounced event deposits that were included and considered in the age-depth model with their calculated ages (Hess et al. 2023b).
Table S3. 14C data of relevant samples from stratigraphic sections at Basta Voe (BV) and Dury Voe (DV) containing tsunami deposits of the Dury Voe event.
The data and sample information were taken from Bondevik et al. (2005) and recalibrated using the most recent IntCal20 calibration curve (Reimer et al. 2020).

Fig. 1 .
Fig. 1.The study area of Loch Flugarth.A. Location of the Shetland Islands at the border between the North Sea and the Norwegian Sea in the North Atlantic.B. Location of the Loch Flugarth in the northernmost part of Mainland.C. The setting of Loch Flugarth and Sand Voe bay, both separated by a ~2-m-high and 70-m-wide barrier.D. Location of Loch Flugarth and position of the sediment cores (FLUG 1-4), and the modern environmental samples M1 and M2.

Fig. 2 .
Fig. 2. Overview of the lithostratigraphy.From left to right: light photograph and computer tomography (CT) scan of FLUG 2; light photograph, CT scan and facies distribution of FLUG 3, including key correlation layers with FLUG 2; high-resolution data of bulk density and magnetic susceptibility; S/Ti and Zr/K ratios (modified from Hess et al. 2023b).The red box indicates the core section shown in Fig. 5.

Fig. 3 .
Fig. 3. A. Factor map showing the results of the hierarchical clustering of principal components (HCPC) with four different clusters and 95% confidence intervals.Numbers indicate the depth of a measurement in cm b.s.(below sediment surface).B. Core FLUG 3 indicating the origin of exemplary samples shown in the HCPC.Colours of the frames refer to colours of the grouping in the HCPC plot (modified after Hess et al. 2023b).PC1 is mainly driven by Br/Ti and S/Ti and, thus, by organic matter.It inversely correlates with heavy mineral content.This may be a function of sediment source from distal (shallow marine to coastal) with low values to increasingly proximal (beach ?catchment ?internal bioproduction).PC2 is mainly driven by Si/Ti and Si/Zr (Fig.S6), reflecting the beach (and dune) sand.

Fig. 4 .
Fig. 4. Boxplots showing the significance of the six most important ratios for each cluster with Bonferroni-corrected p-values.Values on the y-axis indicate how much the four different clusters of the hierarchical clustering of principal components (HCPC) analysis are driven by these individual element ratios.Colours and numbers of clusters are the same as in Fig. 3.

Fig. 5 .
Fig. 5. Core photograph of the basal sand section with sublayers and small peat clasts marked by thin white outlines.The sublayers are shown with mean grain size, sorting, grain-size distribution, magnetic susceptibility (from the present study), bulk density, S/Ti and Zr/K (Hess et al. 2023b).

Fig. 6 .
Fig. 6. rBacon-based age-depth model based on radiocarbon data (blue) and 137 Cs data (green).The posterior age-depth model with the median (red dashed line) and the 2r error range runs from bottom left to top right.The horizontal bars in the figure represent the event deposits >2 mm thick (modified after Hess et al. 2023b).
Recently, Mac Conamhna (2023) reported the possible observation of a tsunami in the earliest written chronicles of the Gaelic world.The chronicle entry from the monastery of Iona in the Outer Hebrides, passed on in three different Medieval annals (Chronicum Scotorum, the Annals of Tigernach and the Annals of Ulster), mentions 'A belch/bursting forth/huge tidal wave/ eruption of the sea in the month of October' (Mac Conamhna 2023).Furthermore, the Annals of Ulster

Fig. 7 .
Fig. 7. Grain-size of the modern beach at Flugarth, and three representative samples showing the range of distributions within one normally graded sublayer of the basal sand layer (cluster 1/facies 4, red lines).Additionally, the distribution of endmember 1 of the endmember modelling analysis in Hess et al. (2023b), representing the thin sand layers from storm overwash above the basal sand is shown for comparison.

Fig. 8 .
Fig. 8. Timeline showing recalibrated 14 C data above and below the tsunami deposits of Dury Voe and Basta Voe (TableS3;Bondevik et al. 2005) in combination with the age of the candidate tsunami deposit from Loch Flugarth based on the age-depth model of the present study (Data S1, Fig.6).Grey curves show the probability distributions, and black or grey vertical lines show the median age.The grey box for Loch Flugarth indicates the maximum and minimum age based on the 2r error.The black box indicates the age interval based on the 1r error.The light blue bar indicates the timing of the assumed seismic tsunami recorded in the Irish annals (Mac Conamhna 2023).The base of the curves refers to the vertical position of the dated sample in reference to the tsunami deposit. 14C dates below the deposits represent maximum ages with an uncertain time gap. 14C dates above the deposits represent minimum ages of the tsunami deposit under the assumption that the dated organic material was embedded into the stratigraphy immediately after death and corresponds to the age of the surrounding deposit.

Fig
Fig. S2. A. The zodiac used for sampling moored on the eastern shore of Loch Flugarth.B. UWITEC gravity corer used for the extraction of FLUG 2-4 (photographs by T. Patel).

Fig. S3 .
Fig.S3.Overview of the lithostratigraphy.From left to right: light photograph and CT scan of FLUG 2; light photograph, CT scan and facies distribution of FLUG 3, including key correlation layers with FLUG 2; highresolution data of Fe/Ti and Sr/Rb ratios, and Sr/Br ratios, respectively (modified fromHess et al. 2023b).The red box indicates the core section shown in Fig.5.The time scale refers to core FLUG 3.

Fig. S4 .
Fig.S4.Photograph and field log of sediment core FLUG 1, taken with a Russian chamber corer at 200-250 cm below lake bottom (b.s.) (Fig.1D, main text).Even though FLUG 2 and FLUG 3 did not reach below the base of the candidate tsunami deposit, FLUG 1 unequivocally shows the lacustrine facies of presumably Holocene age underneath.

Fig. S6 .
Fig. S6.Biplot showing the results of the principal component analysis (PCA).Numbers indicate the depths of a measurement in centimetres b.s.The variables (element ratios) of the PCA are represented by arrows.The length of the arrows indicates the representativeness of each variable in relation to the first two dimensions.

Fig. S7 .
Fig. S7. A. Scree plot showing the percentage of the explained variance per dimension/principal component.B. Quality of representation (cos2) of each variable (proxies) in each dimension.

Fig. S8 .
Fig. S8.Detailed depictions of small rip-up clasts in sublayers 2-4 of the candidate tsunami deposit.A. CT scan in planar view of FLUG 3 at a depth of ~79.5 cm b.s.(sublayer 4) showing an almost round rip-up clast of peat floating in the sandy matrix.B. Several very small rip-up clasts with varying shapes.The largest one is found in the upper part of the normally graded subsequence of sublayer 2. It has an irregular shape and its longest axis has a vertical orientation.Small millimetre-scale peat clasts are distributed in sublayers 3 and 4.

Table 1 .
Facies classification of layers in cores FLUG 2 and 3 based on visual classification and multivariate statistical analyses (principal component analysis (PCA), hierarchical clustering of principal components (HCPC)), and the correlation of both schemes.A brief interpretation of depositional processes is provided for each cluster.