A 1500‐year record of North Atlantic storm flooding from lacustrine sediments, Shetland Islands (UK)

Severe storm flooding poses a major hazard to the coasts of north‐western Europe. However, the long‐term recurrence patterns of extreme coastal flooding and their governing factors are poorly understood. Therefore, high‐resolution sedimentary records of past North Atlantic storm flooding are required. This multi‐proxy study reconstructs storm‐induced overwash processes from coastal lake sediments on the Shetland Islands using grain‐size and geochemical data, and the re‐analysis of historical data. The chronostratigraphy is based on Bayesian age–depth modelling using accelerator mass spectrometry 14C and 137Cs data. A high XRF‐based Si/Ti ratio and the unimodal grain‐size distribution link the sand layers to the beach and thus storm‐induced overwash events. Periods with more frequent storm flooding occurred 980–1050, 1150–1300, 1450–1550, 1820–1900 and 1950–2000 ce, which is largely consistent with a positive North Atlantic Oscillation mode. The Little Ice Age (1400–1850 ce) shows a gap of major sand layers suggesting a southward shift of storm tracks and a seasonal variance with more storm floods in spring and autumn. Warmer phases shifted winter storm tracks towards the north‐east Atlantic, indicating a possible trend for future storm‐track changes and increased storm flooding in the northern North Sea region.


Introduction
Severe storm flooding poses a major hazard to the coasts of north-western Europe with net coastal erosion as well as detrimental effects on both infrastructure and human populations.With rising sea levels as a consequence of climate change, storm-induced coastal flooding is predicted to occur in shorter time intervals in the future (Rahmstorf, 2017;Arias et al., 2021).In north-western Europe, an increase in storm flooding with a poleward shift of storm tracks is expected (Goslin and Clemmensen, 2017;Seneviratne et al., 2021).Simultaneously, global warming reduces the meridional temperature gradient of the high latitudes with the effect that turbulences along the polar front producing the frontal cyclones decrease (Dawson et al., 2007;Sabatier et al., 2012).
To assess whether storm activity has increased in northwestern Europe and to catch high-magnitude and low-probability events, high-resolution records of past North Atlantic storm flooding over timescales longer than the instrumental record are required (Clarke and Rendell, 2009;May et al., 2013;Goslin and Clemmensen, 2017;Moskalewicz et al., 2020;Leszczyńska et al., 2022).These records may be extended by the use of proxy data from sedimentary archives that reflect specific aspects of storm flooding at centennial to millennial timescales, for example from coastal lakes (Liu, 2004;May et al., 2013).
Due to their prominent location in the north-east Atlantic, several long-term records on past storm activity have been established on the British Isles, mostly using aeolian drift sands as storm signals or bromine as sea-spray signal, for example in coastal marshes in England (e.g.Swindles et al., 2018), in machairs (Gilbertson et al., 1999;Dawson et al., 2004b) or in inland ombrotrophic peat bogs (Orme et al., 2015(Orme et al., , 2016;;Stewart et al., 2017;Kylander et al., 2020) of Scotland and Wales.These studies highlight changes in storm frequency during the Medieval Warm Period (MWP, 950-1350 CE) and the Little Ice Age (LIA, 1400-1850 CE) and correlate them with the North Atlantic Oscillation (NAO).The NAO Index reflects the difference between normalised mean winter (December to March) surface air pressure around Iceland and the Azores (Hurrell, 1995).Periods of higher storm activity occurred regardless of whether warmer or transitional phases and a largely positive NAO mode (e.g.Dawson et al., 2002;Orme et al., 2015;Stewart et al., 2017), or colder periods with a negative NAO mode prevailed (e.g.Dawson et al., 2002Dawson et al., , 2004b;;Hansom and Hall, 2009).This association with both phases of the NAO thus provides an ambiguous picture of storm dynamics in the North Atlantic.
The Shetland Islands are located in a crucial geographical position for reconstructing palaeo-storm flooding since cyclones which mostly form near Iceland pass by the Shetlands before reaching the densely populated coasts of mainland Europe (Fig. 1a).Mainland is the largest island of the Shetlands accounting for 65% of the total land area (Bradwell et al., 2019) (Fig. 1b).The only studies of palaeo-storm flooding from northern Mainland are based on dating cliff-top storm deposits, which are quarried by exceptional waves and transported onshore during major storm events.The time of deposition of these events is indirectly derived by dating the peat or sand underlying large boulders, or by lichenometry (Hall et al., 2008;Hansom and Hall, 2009).These chronologies therefore contain considerable uncertainty.Hence, more precise data on the variability of regional storm flooding are needed.
Since they offer precise dating and good preservation of coarser sediments as storm signals (Goslin et al., 2018), peat-dominated coastal environments are among the most suitable sedimentary environments to study storm-induced sedimentary records.The sand in these deposits can be well discriminated from the organic-rich background sedimentation using both sedimentological and geochemical approaches (Osleger et al., 2009).Therefore, to address the knowledge gaps described above, the present study of two sediment cores (spanning ca.1500 years) from the coastal freshwater lake Loch Flugarth in the north of Mainland (Fig. 1c,d) is part of a wider study on the palaeo-storminess of north-western Europe building on existing records and future work.
The aims of this study are to test: (i) if sandy layers from the lake sediment core in the otherwise peaty sediment represent overwash and can be used to reconstruct storm flooding; (ii) if periods of stronger or less intense storm flooding can be deduced from the lacustrine sediment record; and (iii) if temporal variability of past storm flooding can be related to variability of (supra-)regional climate patterns and drivers (e.g. the NAO) to enable future regional predictions.

Geology and geomorphology
The underlying and surrounding bedrock of Loch Flugarth is dominated by schistose feldspathic quartzites and thickbanded hornblendic and granulitic gneisses of the Sand Voe group.The gneisses contain up to 50% plagioclase (K-feldspar) and the quartzites are partly impured by muscovite and biotite (Pringle, 1970).
The Shetland Islands were glaciated during the Pleistocene (Bradwell et al., 2019;Hall et al., 2021).Postglacially, relative sea level (RSL) rose more or less continuously, but has slowed down since the mid-Holocene, with RSL rise probably in the range of only a few centimetres over the last 1500 years (Bondevik et al., 2005;Dawson et al., 2006).Postglacial emerged beaches or raised glacial-marine sediments as an indication of a mid-Holocene RSL highstand are absent in Shetland (Kelley et al., 2018;Hall et al., 2021).
The coastline of the Shetland Islands varies greatly, ranging from cliffs of up to 300 m height to peat-floored beaches (Flinn, 1974).Beaches and dunes are rare due to low sediment availability and occur only in very sheltered locations (Hall et al., 2021), such as at the study site of Loch Flugarth and Sand Voe bay.Peat is the major component of Shetland's landscape today, covering ca.50% of the entire archipelago (Hall et al., 2021), mostly exceeding a thickness of 100 cm (Dry and Robertson, 1982).Its formation presumably began during the Windermere Interstadial from 14.7 to 12.9 ka BP (McIlvenny et al., 2013).

Climate
The Shetland climate is oceanic temperate, humid, with long mild winters and short cool summers (Köppen-Geiger climate class Cfc; Beck et al., 2018).Annual mean temperature in the capital city Lerwick is 7.4 °C, and ranges from 3.6 °C in February to 12.4 °C in July and August (related to mean monthly temperature, period 1930-2021;Met Office, 2021).Due to the warming influence of the North Atlantic and the Shelf Edge Current to the west, the climate is warmer than in other areas of comparable latitude (Bigelow et al., 2005;Kelley et al., 2018).
Precipitation is mostly associated with the passage of Atlantic depressions, which are prominent during autumn and winter (Met Office, 2016).Autumn and winter are slightly wetter with monthly mean precipitation of 146 mm in November compared to 53 mm in May.Annual mean precipitation is 1254 mm (Met Office, 2021).
Shetland is one of the windiest places of the British Isles with an annual mean wind speed of 14.6 knots (kn) and 58 gale days per year (Irvine, 1968;Bennett et al., 1992;Hall et al., 2021).Mean wind speed is higher during winter (Met Office, 2016) with the highest mean wind speeds recorded in January (18.2 kn) and the lowest in July (10.9 kn) (Irvine, 1968).In winter, gales occur on average in two days out of five and storms can arrive from varying directions, but mostly originate from the north, west and south-west (Birnie et al., 1993;Met Office, 2016).Atlantic depressions originating in the west reach the British Isles and generate winds from southern and south-western directions.Subsequent to passing the British Isles, the wind direction switches to north or north-west, which typically generates the strongest gusts (Lamb, 1991;Met Office, 2016).In summer, Atlantic storms take more northerly tracks and are usually weaker (Irvine, 1968).
There are three different storm track types (Fig. 1a): (i) the fast-travelling Jutland-type cyclones with high wind speeds originating over Newfoundland and passing Jutland; (ii) the Scandinavia-type cyclones created over Greenland and Iceland crossing Scandinavia north of 60°N; and, (iii) the Skagerrak type with tracks running between types i and ii.The Scandinavia-type cyclones have long dwell times and sometimes appear almost stationary, and therefore often coincide with high tides (Petersen and Rohde, 1991;Wöffler, 2016).
Storm surges on the Shetland Islands are particularly driven by the Scandinavia and the Skagerrak types.

Loch Flugarth and Sand Voe bay
Loch Flugarth (60°36′N, 1°20′W) is a small (0.16 km²) lake located in north-western Mainland (Fig. 1b).The lake level is about 2 m asl (above mean sea level) and separated from Sand Voe bay by a ca.70-m-wide sand-and-gravel barrier (Fig. 1c).The barrier has a flat beach and is topped by a dune reaching ca. 2 m above mean high tide level.Sand Voe is a relatively shallow bay opening towards the north-west and is framed by cliffs higher than 10 m (Flinn, 1974).A cemetery and a road are located between the sand beach and the lake indicating low mobility of the barrier system in recent times.
Loch Flugarth has a flat bottom, a maximum water depth of ca.2.4 m (Murray and Pullar, 1909) and a low alkalinity due to the prevailing grassland, heather and bog surrounding it (UK Centre for Ecology & Hydrology, 2004).It is characterised by the absence of thermal stratification (Bennett et al., 1992).Two small inflows from (south-)western directions are the only tributaries into the lake.There is an outflow between the beach and the lake.During high tides the beach berm is entirely submerged by seawater and the foredunes show erosional scarps from recent winter storms (Supporting Information Fig. S1).This indicates that low-frequency, high-intensity storms may create overwash and introduce seawater and sand from the beach and dunes into the lake as described, for example, in Liu (2004) and Chaumillon et al. (2017).

Potential limitations of the study site for reconstructing storm overwash
Accurate reconstruction of past storm overwash is based on the prerequisite that the sedimentary archive did not change significantly.Available satellite imagery and historical maps do not show any notable changes during the last 200 years.RSL rise would increase the archive's sensitivity by shifting the shoreline inland (Sabatier et al., 2012) and lowering the threshold of overwash (Liu, 2004;Biguenet et al., 2021).However, Bondevik et al. (2005) have shown that significant RSL changes over the last 1500 years are unlikely.This is also supported by a combination of recent measurements and geological proxy data, including continuous Global Positioning System measurements, glacial isostatic adjustment models and tide-gauge data from Shetland (Wahl et al., 2013).
Non-varved polymictic lakes with a simple bathymetry such as Loch Flugarth are less prone to internal lake processes such as hyperpycnal flows or interflows and are therefore very suitable as climate archives (Schillereff et al., 2014).However, water-column mixing and water-level fluctuations are more likely in shallow lakes, which may affect continuous sedimentation.Wind-generated waves on the lake surface may generate turbulence in the water column and promote resuspension of bottom sediment (Schillereff et al., 2014).However, the lake has a very small fetch and wind waves are expected to have only a minimal impact on the coring site.In addition, post-depositional erosion as well as disturbances from erosion at the short subaquatic slopes can be excluded, as both cores FLUG 2 and 3 were taken in the deepest parts.

Fieldwork
Four short sediment cores were taken from the northern part of Loch Flugarth at a water depth of ca.2.0 m (Figs.1d and S1) in March 2018 using a zodiac and a UWITEC gravity corer, equipped with transparent PVC liners with an inner diameter of 60 mm (Fig. S2b).Three short sediment cores (FLUG 2-4) were taken from the lake bottom down to depths of up to 91.7 cm (Fig. 1d).Core FLUG 1 was taken using a Russian chamber corer mounted on a manual Edelman coring system and opened at a depth of 200 cm below the lake bottom in order to sample at a depth of 200-250 cm below the lake bottom.Additionally, two surface samples of modern sedimentary environments were taken from the shallow lake bottom close to the barrier (FLUG-M1) and from the uppermost subtidal zone of the marine embayment (FLUG-M2).In this study, we focus on sediment cores FLUG 2 and 3 with a length of 77.0 and 91.7 cm, respectively, and the two modern surface samples (FLUG-M1 and -M2).

Physical and textural sediment analysis
Sediment cores were cut and split in the laboratory of the Geological Survey of Belgium.The lithostratigraphy was described in terms of grain size, colour and organic macroresidues.In total, 91 samples from FLUG 3 were taken with an average step size of 1.0 cm.At some depths, the step size was smaller to stay within lithostratigraphic boundaries.The samples were taken in 2.0-cm steps from 79 cm below the sediment surface (b.s.).

Geotek core logging and computer tomography
Core FLUG 3 was analysed using a Geotek multi-sensor core logger at the Renard Centre of Marine Geology, Ghent University, Belgium.The core logger was equipped with a Geoscan IV line scan camera for high-resolution core photographs, a γ-ray attenuation densometer to measure bulk density and a spectrophotometer (Konica Minolta CM-700d) to measure the spectrophotometric reflectance and the CIE colour spectrum (Geotek, 2021).All parameters were measured with a spatial resolution of 2.0 mm.X-ray computer tomography (CT) scans of FLUG 2 and FLUG 3 were produced by a Siemens Somatom Definition Flash Medical X-ray CT scanner at Ghent University Hospital.The reconstructed image stacks have voxel sizes of 0.15 × 0.15 × 0.3 mm.Stacked images of the cores were analysed using the Fiji image processing package (Schindelin et al., 2012).The scanner creates a three-dimensional (3D) space to analyse the internal structure of the sediment core.

Grain-size analysis by laser particle diffraction
Grain-size analysis was carried out in the Laboratory for Geomorphology and Geoecology, Heidelberg University, Germany, using a Fritsch Analysette 22 NeXT Nano device with a measuring range of 0.01-3800 μm (Fritsch GmbH, 2020).The pre-treatment of samples included drying at 105 °C and manual grinding to carefully break up aggregates.Organic matter was dissolved using H 2 O 2 (30%).Aggregates were further dispersed using a distilled water solution with 0.1 mol Na 2 P 2 O 7 (44.6 g l -1 ) in an overhead shaker.Triplicate measurements of grain-size distribution (GSD) were averaged and univariate statistical measurements were calculated after Folk and Ward (1957) using GRADISTAT (v.9.1) (Blott and Pye, 2001).To statistically group entire GSDs and correlate them with facies types, end-member modelling analysis (EMMA) was applied using the R package EMMAgeo (Dietze and Dietze, 2019;Dietze et al., 2012).

Geochemical sediment analysis
Total organic and inorganic carbon Total organic carbon (TOC) and total inorganic carbon (TIC) were analysed using an Elementar soliTOC cube in the Laboratory for Geomorphology and Geoecology, Heidelberg University, Germany.This analysis is based on dry combustion with a temperature ramp procedure and CO 2 detection with a detection limit of <15 μg (Elementar, 2020).

X-ray fluorescence (XRF) spectroscopy
Non-destructive XRF spectroscopy was performed using an Avaatech (GEN-4) XRF core scanner at the Institute of Earth Sciences, Heidelberg University, Germany.The core scanner was equipped with an OXFORD 'Neptune 5200' series 100 W X-ray source with a rhodium anode.The slit size was adjusted to 5.0 mm downcore and 10.0 mm cross-core.Measurements were carried out at two energy levels: specifically, at 10 kV with a current of 500 μA without a filter and a dwell time of 10 s and at 30 kV beam with a 1500 μA current, a Pd-thick filter and dwell time of 10 s.Data processing was performed using bAxil software (Brightspec, 2015).To remove potential elemental bias due to topographic variations of the core surface, such as water content or porosity, the elements were normalised by dividing the intensity counts of each element by the XRF total during each run excluding coherent and incoherent scattering (Kern et al., 2019).
Interpretation of the geochemical analysis is based on a principal component analysis (PCA), which was carried out on 13 centred log-ratio-transformed elements (Al, Br, Ca, Fe, K, Mn, Mo, P, Rb, S, Si, Sr, Zr).The hierarchical clustering of principal components (HCPC) was applied based on the first five PCs of the PCA to discriminate the most meaningful variables (i.e.elements) and to group similar sediment layers with overlapping characteristics from the variables.The PCA and the HCPC were performed in R studio following Kassambara (2017) using the R packages FactoMineR (v.2.4) for computing and factoextra (v.1.0.7) for visualisation of the results.

Chronology
A combination of 14 C and 137 Cs dating was used and integrated into a Bayesian age-depth model.The 14 C samples of both cores, FLUG 2 and 3, were taken at six different depths and analysed by accelerator mass spectrometry (AMS) in the Poznań Radiocarbon Laboratory, Poland (Goslar et al., 2004).Five 14 C dates were generated from terrestrial plant fragments and one from bulk peat, as plant macrofossils were extremely rare.All 14 C samples were taken from FLUG 3 except for one sample from FLUG 2 (Table 1).The 14 C data were calibrated using the IntCal20 curve (Reimer et al., 2020).A lake-reservoir effect based on older peat introduced by surface runoff or by old groundwater in similar environments (Björck and Wohlfarth, 2001) was circumvented by dating macro-plant fragments in all cases except sample FLUG 3-30-31 (Table 1).
For 137 Cs dating, the upper 15 cm of one half of core FLUG 2 was sampled at 1-cm resolution.Material was freeze-dried and homogenised at the Institute of Earth Sciences and in the Laboratory for Geomorphology and Geoecology, Heidelberg University, Germany.The 137 Cs depth distribution of the FLUG 2 sediment core was measured at Ravensburg-Weingarten University and modelled by three 137 Cs peaks (Chernobyl: 1986;nuclear weapons tests: 1963 and1959) using a model of Smith (Aquascope) (Smith et al., 2005) (Fig. 2; Table S1).Each peak decreases with a sum of three exponential functions, each for wash-off, increasing fixation and equilibrium processes.The onset of radioactive fallout by fusion bomb tests in 1952 is also visible in the 137 Cs depth distribution and was roughly approximated in the model by a constant value between 1952 and 1963.Considering diffusion and turbation, the sum of the nine exponential functions and the constant were convolved by a Gaussian curve.Parameters of the model were optimised with a least-squares fit procedure.
Bayesian age-depth modelling was carried out using rBacon (Blaauw and Christen, 2011).It was based on the calibrated 14 C and 137 Cs data and considered the influence of sedimentation conditions, for example changing sedimentation rates or events.A depth file of 0.2-cm steps was created and implemented into the model to better correlate the resulting years with the Geotek and XRF data.A thickness of 1 cm was set resulting in an age-depth model with 70 sections to obtain a model as smooth as possible with sufficiently high resolution (Blaauw and Christen, 2011).

Evaluation of historical records
In an attempt to correlate individual sand layers of the sediment cores with specific overwash events such as major storms, pre-instrumental data were exploited from historical documentation from Lamb (1991).Lamb (1991) describes and categorises over 166 storms for the entire North Sea region using scientific literature, local archives, diaries and newspapers since 1500 CE.The author establishes a storm severity index based on the rating of the following parameters: 'greatest wind speed', 'greatest area covered at any stage by winds causing widespread damage' and 'total duration of the occurrence of damaging winds'.The storm severity index ranges from Class VI (low severity) to Class I (very high severity).Four conditions were established to consider an historical storm to be relevant for the sedimentary archive of Loch Flugarth: i. minimum Class V, ii.minimum Beaufort Scale 10 (if available), iii.geographical description with tags: 'Shetland', 'Northern North Sea', 'Northern Scotland', 'Northern part of the British Isles' or 'Eastern North Atlantic', and iv.western to northern wind direction, due to the exposure of Sand Voe's entrance.
All considered storms were chronologically compiled in an Excel file with all relevant information such as date, time span,  duration, wind direction and speed, severity (after Lamb, 1991), location, description and reference (Hess et al., 2023).All dates were expressed in the Gregorian calendar and historical dates before October 1582 CE were stated in the new style.

Chronostratigraphy
Age-depth modelling revealed a chronostratigraphy of about 1400-1500 years, dating back to 426-787 cal a CE.The core resolution across the entire core was 20 a cm −1 , whereas the uppermost 15 cm dated by 137 Cs revealed a higher core resolution of about 10 a cm −1 (Fig. 2).Further details on the results of age model construction can be found in Supplementary Note 1 and Table S2.

Lithostratigraphy
Cores FLUG 2 and 3 showed lacustrine sediments being visually separated into four different facies.The two cores were taken ca.20 m apart.Their stratigraphic correlation was excellent (Fig. 3).All lithostratigraphic descriptions refer to depths identified in FLUG 3. Facies 1 comprised muddy, organic-rich, dark brown sediment.At 37.5-43.0and 73.0-77.0cm b.s., respectively, these organic-rich deposits were much brighter in colour and defined as Facies 2. Facies 3 comprised muddy, organic material mixed with sand grains, dominant in the lower part of the core at 57.0-73.0cm b.s.Facies 4 was represented by bright sand layers of varying thickness from a few millimetres up to several centimetres.At the base of the core, Facies 4 was ca.13.0 cm thick (Fig. 3).These basal sands did not represent the bottom of the lacustrine sequence as in FLUG 1 (200-250 cm b.s.; Fig. S3) dark, organic-rich background deposits of Facies 1 still continued.
Two pieces of wood were identified: one of 4.0 cm length at 11.0-12.0cm b.s. in FLUG 2 (Fig. 3a) and another of 3.5 cm length at 12.0-13.5cm b.s. in FLUG 3 (Fig. 3b).A rounded fine pebble (>0.5 cm) and a flat, subangular, grey pebble was found at a depth of 18 and 30.5 cm b.s., respectively (Fig. 3c,d).All pebbles in the core were found at upper or lower contacts of sand layers and appear to originate from the modern beach (FLUG-M2).
The density data of FLUG 3 (Fig. 4) reflected facies changes well, with the lowest densities (<1.0 g cm -3 ) found in the upper 3 cm and the highest densities in the lighter sandy layers (peak densities of about 2.0 g cm -3 ).The CIE colour space data reflected similar trends, such as *b axis (yellow colour space, proxy for light sand).Most measurements ranged between -1.0 and +2.5.Only at depths with purest sand layers, such as 35.9-36.3,75.9-76.5and, in particular, 81.7-91.7 cm b.s., were values >3.0.

Sedimentological characteristics
The grain-size analysis included particles from 0.01 μm (clay) to 700 μm (coarse sand) (Fig. 4).The GSD ranged from unimodal and well-sorted to trimodal and very poorly sorted.Sand content >90% was detected at depths of 4. In general, Facies 4 had a unimodal, poorly to moderately well-sorted GSD; it was better sorted than Facies 1. Modern samples from the lake margin close to the sandy barrier (FLUG-M1) and the uppermost subtidal zone (FLUG-M2) had very similar grain-size characteristics (Fig. S4).Both GSDs were unimodal and moderately well-sorted fine sand with sand content of >96% and a mean value of 218 μm.
End-member modelling of GSDs revealed three distinct endmembers (EM) (Fig. 5).EM 1 depicted a unimodal, moderately well-sorted symmetrical curve with a mode in fine sand.It was similar to some of the samples from the basal sand unit, but seemed to represent a mixture between different sources and processes.EM 3 was a very poorly sorted, trimodal curve with peaks at about 5, 20 and 200 μm, categorised as a fine sandy coarse silt.It had a large overlap with Facies 1-3, in particular samples from 34.0-35.0 and 57.5-58.5 cm b.s.Both surface samples from the seaward barrier and the uppermost subtidal zone resembled EM 1 (Fig. 5).

Geochemical proxies
PCA was first conducted on the elemental concentration with the first two PCs explaining 47.2% (PC1) and 14.0% (PC2) of the variance of the whole dataset (Fig. S5).PC1 was dominated by the variance of Al, Br, K, Rb, S, Si and Zr, and PC2 was dominated by Ti.
XRF results are typically presented in ratios with conservative elements (e.g.Ti, Fe, Si) as a denominator.The ratios were selected based on the results of the PCA of the elements (Fig. S5).Ti was used in particular, as it is not involved in biogenic processes and is not affected by either redox conditions or diagenetic overprinting (Chagué, 2020).Based on the eigenvalues, five dimensions and four clusters were collectively considered for the HCPC.In the HCPC (Fig. 6a), Facies 4 was geochemically divided into basal sands (blue cluster 1), represented by Zr/K and Sr/Br, and upper thin sand layers (yellow cluster 2), mostly driven by Si/Ti and K/Ti (Fig. S6).The remaining darker Facies 1 sediments were separated into a pure, fine-grained, organic-rich group (red cluster 4), for example at 20.0-30.0 cm b.s., with high concentrations of Ca, S, Br and Mn, and a mixed group (grey cluster 3), mostly comprising samples from Facies 2 and 3 (Fig. 6b).
The most representative ratios with high values for each cluster were Zr/K (cluster 1), Si/Ti (cluster 2) and S/Ti (cluster 4) (Fig. 7).The TOC content of FLUG 3 varied between nearly 0% at the bottom of the core and up to >35.0% at a depth of 29.0 cm (Fig. 7).The measured TOC was higher in Facies 1 layers with an average of about 20.0%.The lowest values occurred in homogenous sand layers (e.g. in L41-42).TOC of FLUG-M1 and FLUG-M2 was nearly 0%, too.TIC values were close to 0% in the entire core.

The historical record
The evaluation of storm flooding from historical data revealed that a total of 77 storms affected northern Scotland and the northern North Sea originating from a western to northern direction (Lamb, 1991).Based on the total number of storms, a high storm frequency was recorded in the second half of the 20th century, particularly during the 1980s, and in the 19th century, with seven and eight severe storms during the 1880s and the 1820s, respectively (Fig. 8).For the first half of the 20th and the 18th century, a comparably low number of storms was documented.In contrast, the 1690s-1710s were known for their intense storms with recurring sandstorms along the Scottish coasts (Lamb, 1991;Bampton et al., 2017).A longer period of high storm frequency can be observed between the 1540s and 1610s with a calmer phase in the 1560s (Lamb, 1991).
Most of the storms occurred in winter (Dec, Jan, Feb), particularly those of the 20th century (Fig. 9).The 19th century was characterised by several heavy summer storms (Jul, Aug, Sep) and severe winter flooding due to higher rainfall.Most notable is the high percentage of storms in spring (Mar, Apr, May) during the 18th century.During the 16th and 17th centuries, storm activity was high during late summer and autumn (Aug-Nov).Thus, there is an apparent shift from spring and autumn storms during the 16th-19th centuries towards a dominance of winter storms in recent times (20th century) (Lamb, 1991).

Sedimentary facies interpretation
The EMMA revealed tri-to bimodal GSD curves for EM3 reflecting Facies 1.Despite thorough continuous H 2 O 2 treatment over several days, small amounts of organic matter may have remained in the sediment, which could explain the third mode in the fine silt category at about 4 μm.The high organic-matter content originated either from autochthonous algal productivity or terrestrial run-off from peaty soils of the lakeshore and catchment (Biguenet et al., 2021).The latter may indicate heavier precipitation and/or greater snow melting (McIlvenny et al., 2013).Otherwise, the fine-grained fraction could have entered as fine dust by wind.
Facies 1 was rich in Br, Ca and S.Although Br has been used as marine indicator in other studies from Scotland (Orme et al., 2015;Stewart et al., 2017), Br, Ca and S are also indicators for fine-grained organic matter (Chagué, 2020;Biguenet et al., 2021)  biogenic production in the lake or originate from marine shell fragments introduced from the beach and shallow subtidal zone (Kelley et al., 2018).However, neither shell remains nor any other carbonates were found in the modern intertidal environment (samples FLUG-M1 and FLUG-M2), nor do carbonates form in an acidic, peaty environment such as Loch Flugarth.Kylander et al. (2020) found high Ca concentrations in peat-bog samples in southern Scotland and related them to the underlying metagabbro.The surrounding bedrock of Loch Flugarth is dominated by quartzites and gneisses (Pringle, 1970), which are not known for high Ca concentrations.A calcareous fixed dune north-east of Loch Flugarth on top of the sand-and-gravel barrier (Dargie, 1998) (Fig. 1c) could be a source for Ca introduced by wind or surface runoff during heavy precipitation.Facies 2 was very silty whereas Facies 3 was poorly sorted and had no clear elemental signal, probably representing a mixture of different sediment types.Both Facies 2 and 3 could be assigned to the mixed cluster 3 in the HCPC.In contrast to Facies 1-3, Facies 4 was better sorted and showed a unimodal GSD, best explained by EM1 in the EMMA peaking in the sand category.Additionally, EM1 was nearly identical to the GSD type for the modern samples FLUG-M1 and FLUG-M2, pointing to an unequivocal sand source at the beach and the shallowest part of Sand Voe bay.The lowermost sand layers L41 and L42 deviated from EM1; their GSD is explained to a larger percentage by EM2.This discrepancy and the formation of L41 and L42 are discussed in detail in Engel et al. (2023).
Facies 4 was represented by cluster 2 in the HCPC, dominated by the ratios Si/Ti and K/Ti.In Facies 4, K and Si were the principal elements in quartz and K-feldspars, which are the typical components of the underlying and surrounding bedrock of Loch Flugarth (Pringle, 1970) and commonly high in beach sand (e.g.Pouzet and Maanan, 2020).Si/Ti and K/Ti have been used to differentiate between coastal beaches and river-sand populations (Goslin et al., 2018) and as storm indicators in peat bogs (Kylander et al., 2020).The overall high Fe content in all sediments of the core could be attributed to the humus-iron podzols surrounding Loch Flugarth (Dry and Robertson, 1982).
Both aeolian transport and marine overwash may be responsible for the deposition of thin sand layers in the lake since both processes generate well-sorted GSD curves (Dawson et al., 2004b;Morton et al., 2007;Moskalewicz et al., 2020).In a study on coastal sand movements in the Scottish Outer Hebrides, well-sorted sediments with unimodal curves and a mean grain size of about 200 μm, very similar to Facies 4, were linked to aeolian processes, because marine indicators such as shells or other marine microfossils were absent (Dawson et al., 2004b).In addition, North Atlantic cyclones are associated with high rainfall although wet sand may be also very efficiently transported by wind, for example by saltation (Rotnicka, 2011).However, the sharp, possibly erosional lower contact of some sand layers (Moskalewicz et al., 2020), for example L3, L5,or L12, and the pebbles and wood fragments directly found above or below a sand layer rather indicate storm overwash.
A period of higher storm flooding from 1150 to 1300 cal a CE was in agreement with a proxy record from northern Scotland with a strong presence of Br for the MWP, which was interpreted as a wetter and stormier period (Stewart et al., 2017).Orme et al. (2016) provided evidence for higher storm activity in two ombrotrophic peat bogs on the Outer Hebrides dated to ca. 1150 and 1550 cal a CE, which might correlate with L30-33 and L18-21 in the Flugarth record, respectively.The 14th, 17th and 18th centuries were predominantly calmer periods reflected by fewer and thinner sand layers.The conspicuous sediment facies change at ca. 40 cm b.s.coincided with the period 1400-1420 cal a CE, which has been described as a strong North Atlantic climate change by Dawson et al. (2007).
Several studies have attributed individual sand layers to specific historical storm events during the last 300 years (Sabatier et al., 2008;Moskalewicz et al., 2020;Biguenet et al., 2021).The correlation with historical data, however, was challenging due to the high number of severe storms around the Shetland Islands.For instance, more than 21 severe storms occurred between 1950 and 1990 (Lamb, 1991), but only seven sand layers were identified in the sediment cores for the same period.This discrepancy is quite common as only the most intense storms tend to be preserved in sediment records (e.g.Liu, 2004;Dezileau et al., 2011;May et al., 2013;Moskalewicz et al., 2020).Furthermore, two storms of identical intensity and characteristics may not necessarily cause the same marine flooding and sediment pattern.Changing boundary conditions such as tides, fetch, season, etc., additionally influence the occurrence of overwash (Liu and Fearn, 2000;Liu, 2004;Moskalewicz et al., 2020).Also, the stage of coastal landform development has to be considered (Leszczyńska et al., 2022).Any estimation of the number of storm events from individual sand layers should therefore be treated with caution (Biguenet et al., 2021).Nevertheless, the high abundance of sand layers within the Flugarth records implied a good preservation potential of the record especially when compared to other studies (Sabatier et al., 2008;Pouzet and Maanan, 2020;Biguenet et al., 2021).
Layer L1 most probably corresponded to a storm in the late 1980s or 1990s, possibly the severe storms in 1992 or 1993, where boulders at Eshaness on the north-west coast of Mainland close to the study site were shifted on top of the >40-m-high cliffs by waves (Hall et al., 2006;Hansom and Hall, 2009).Layer L2 might represent the severe storm of December 1986, while L3 could be correlated with a storm either in 1983 or in 1981 as the 137 Cs peak of 1986 lay in between L2 and L3 (Fig. 2b).The historical North Sea storm in 1953, which caused a 0.5-3.0-cmthicksandy layer in coastal peats in eastern England (Swindles et al., 2018), was also described as severe for the Shetland Islands and was probably represented by one of the sand layers L5-7.
The pronounced sand layers L11-12 might reflect the severe storms in the middle of the 19th century, for example the storms in 1839 or 1855 (Lamb, 1991).The relatively thick layer L12 at 16 cm b.s.dated to 1735-1901 (median: 1831) cal a CE might have been caused by a massive storm reconstructed for North Scotland (McIlvenny et al., 2013).Although large sand movements with accumulation of several decimetres were reported in the 1690s in Scotland and the southernmost part of Shetland, for example the Udal storm or the Culbin Sands Disaster (Lamb, 1991;Orme et al., 2016;Bampton et al.,

Atmospheric and oceanic drivers of storm overwash
Both the sedimentological and the geochemical results (Fig. 10a,b) have shown that periods of more frequent storm overwash in the Loch Flugarth sediment core occurred during the MWP until ca.1300 and 1450-1550 cal a CE.After a somewhat calmer period during the coldest stages of the LIA, storm overwash intensified after the second half of the 19th century with an intermittent period of lower storm flooding during the first half of the 20th century, which is also in line with the historical data (Fig. 10c).
North Atlantic water inflow based on coccolith populations (Giraudeau et al., 2010) was considered to have a major  (after Lamb, 1991).Storm descriptions missing information about the month of the storm were not included, which resulted in data gaps especially in the 16th and 17th centuries.(Orme et al., 2016).(e) Cliff-top storm deposits on the Shetland Islands (Hansom and Hall, 2009).(f) Winter (DJF) NAO reconstruction (Trouet et al., 2009).(g) Smoothed Greenland sea salt influences (ssNa+) of core GISP 2 (Meeker and Mayewski, 2002).(h) Stalagmite proxy from a Scottish cave (Baker et al., 2015).(i) Reconstruction of Atlantic water inflow in the Norwegian Sea derived from concentrations of the coccolith Gephyrocapsa muellerae (Giraudeau et al., 2010).[Color figure can be viewed at wileyonlinelibrary.com] climatic impact on the Loch Flugarth study area due to the proximity of the Shelf Edge Current to the north-west coast of Shetland.However, except for the last 150 years, there is no clear correlation between increased warm water inflow from the tropics and higher storm flooding in Flugarth (Fig. 10i).
A proxy of annually laminated stalagmites from a Scottish cave based on lamina counting and U-Th dating (Fig. 10h) reflected the transition from the warm MWP into the colder LIA around 1400 CE (Baker et al., 2015) and supported the persistent positive NAO mode during the MWP (ca.1000-1400 CE), presented by Trouet et al. (2009) (Fig. 10f).Storm overwash in Loch Flugarth was also high during the MWP.However, there were calmer storm flooding periods in between from ca. 1300 cal a CE until ca.1450 cal a CE which do not fit with a positive NAO mode.During the LIA, the NAO mode was mostly negative with short periods of a positive NAO (Fig. 10f,h), which could be explained by high-pressure systems that formed over the north-east Atlantic due to colder surface waters weakening the Icelandic deep pressure system.This was attributed to stationary high-pressure systems which would hinder the development of cyclones in the North Atlantic and thus reduce storm activity during winter (Orme et al., 2016).
This discrepancy in records may be resolved by considering seasonality patterns (Trouet et al., 2012).Since the NAO is driven by winter weather patterns, the higher storm flooding during the negative NAO mode (i.e.LIA) may be explained by severe cyclones predominantly occurring in spring and autumn (Lamb, 1991;Trouet et al., 2012).Additional evidence for this hypothesis is provided by Wheeler et al. (2010), who inferred severe storms during spring and autumn between 1685 and 1699 from Royal Navy logbooks.This is in line with the increased storm frequency observed in the historical evaluation of storm frequency per century in spring and autumn during the 17th to 19th centuriesy (Lamb, 1991) (Fig. 9).However, this seasonal shift does not explain why these storms were not well recorded in the Loch Flugarth record.
One hypothesis could be a different prevailing wind direction.Based on model simulations, Bampton et al. (2017) stated that present wind fields are similar to the wind conditions of storms observed during the LIA.However, historical documents indicated more easterly and north-easterly winds during the LIA (Lamb, 1991), which would probably produce less overwash at Flugarth due to the northern-western to northern opening of Sand Voe bay.In addition, the lake sediment layers with high Ca content were dated to the LIA.North-easterly winds during the LIA could have transported Ca from the calcareous dune immediately north-east of Loch Flugarth (Dargie, 1998) (Fig. 1c), which under a colder and drier climatic regime (Stewart et al., 2017) might have been less vegetated, into the lake.The high peat content in the sediment core during the late LIA was possibly introduced by increased surface discharge associated with precipitation maxima in summer or by shoreline peat erosion during the 18th and 19th centuries (Dawson et al., 2004a).
Another hypothesis could be that the southward extension of the Arctic sea ice during the LIA may have cooled the air and the sea surface in the north-east Atlantic (Dawson et al., 2002(Dawson et al., , 2010)).Consequently, the storm tracks were shifted further south (Dawson et al., 2002(Dawson et al., , 2010;;Giraudeau et al., 2010;Orme et al., 2015;Stewart et al., 2017).Biguenet et al. (2021) emphasised that only storm tracks which ran very close on the seaward side of the coast could be captured in their sedimentary archive of a small coastal lagoon.Thus, even a slight shifting of the storm tracks could determine whether a storm of equal intensity was preserved in the lake record.Since the site only recorded storms with tracks north of 60°N, Jutland-type cyclones moving from south-west to north-east were therefore unlikely to be captured by the Flugarth record (Fig. 1a).The hypothesis of a larger-scale southward shift of storm tracks is supported by comparisons of studies in the Mediterranean, which found a clear link between cold phases of the Holocene (e.g.around 1550-1990) and higher storm frequency, whereas storm frequency was low during the MWP (Sabatier et al., 2008(Sabatier et al., , 2012;;Costas et al., 2012).
Cliff-top storm deposits (Fig. 10e) were found at 15-60 m asl in Shetland, and dated to 700-1050, 1300-1900 and post-1950 cal a CE (Gilbertson et al., 1999;Hansom and Hall, 2009).Due to the proximity of the cliff-top storm deposit site to Loch Flugarth and a similar orientation of the cliffs, the high storm flooding between 1300 and 1900 should be also visible in the Flugarth record.However, Hansom and Hall (2009) acknowledged that the dating precision of the cliff-top storm deposit is relatively low due to indirect dating of adjacent sand or peat underneath the boulders.Thus, the cliff-top storm deposit record during the LIA might also relate to a low number of high-intensity events instead of a generally higher storm frequency (Trouet et al., 2012;Orme et al., 2016;Stewart et al., 2017).The higher intensity might be explained by the generally steepened thermal gradient between the cold polar waters and the warm subtropics during colder periods (Dawson et al., 2007;Sabatier et al., 2012).Thus, it could be inferred that the storm frequency was relatively low during the LIA (Orme et al., 2016), resulting in a smaller number of events coinciding with high tides and a lower probability of overwash at Flugarth.
Further arguments for a moderate storm overwash during the LIA are that aeolian sand drift does not only depend on wind strength, but also on sand availability, vegetation cover, anthropogenic impact and drier conditions (Bampton et al., 2017).The coldest stages of the LIA were associated with a reduced winter storm frequency over the North Atlantic due to the expansion of the polar anticyclone, providing a negative NAO mode and stable conditions during winter (Dawson et al., 2004b).These highpressure systems could have created drier conditions during the LIA (Stewart et al., 2017) and, hence, reduced vegetation cover.Together with anthropogenic landscape changes (i.e.overgrazing), sand dunes might have been less stable (Bampton et al., 2017).The recent decline in aeolian sand transport observed in Scotland over the past 200 years (Fig. 10d) could be attributed to anthropogenic adaptation activities on dune systems (Orme et al., 2016).
Our study is a single-site study reflecting storm flooding under local conditions.Given the ambiguities of single records pointed out in Wallace et al. (2021), further regional studies from the northern North Sea and in particular the Shetland Islands are required to better filter out local effects of the Flugarth site and complement the pattern of north-east Atlantic palaeo-storm dynamics.

Conclusions
This study has shown that the two lake sediment cores from Loch Flugarth in the north-western part of Mainland, Shetland Islands, recorded the local North Atlantic storm flooding over the last ca.1500 years.This multi-proxy approach using sedimentological, geochemical, chronological and historical data revealed 42 sand layers originating from the beach and the very shallow subtidal zone.The marine origin of the sand layers could be concluded from the grain size characteristics and the XRF-based Si/Ti and K/Ti element ratios discriminating these sand layers from the other sediment types.
The data clearly showed periods of higher storm flooding around 980-1050, 1150-1300, 1450-1550, 1820-1900 and 1950-2000 cal a CE.However, a clear correlation of individual sediment layers to single historical storm events was difficult due to the high number of historically documented storms.
Most periods of higher storm flooding correlated partly with positive NAO phases and enhanced North Atlantic water inflow into the Norwegian Sea.It can be inferred that in the phases when no correlation takes place, other, local factors, such as wind direction and vegetation cover, may have influenced the frequency and intensity of storm overwash at the study site.The number and thickness of storm layers dated to the MWP, particularly 980-1050 and 1150-1300 cal a CE, as well as during the second half of the 20th century, suggested similar storm conditions.The lake sediments reflected storm tracks running north of the Shetland Islands and storms presumably originating from a north-westerly to northerly direction due to the orientation of Sand Voe bay.Additionally, a shift of storm tracks towards the north-east Atlantic during warmer periods was inferred.
The LIA was dominated by a different wind regime and storm tracks probably shifted further south, resulting in higher storm flooding in southern parts of Europe.Storm flooding shifted both spatially and temporally with more storms in spring and autumn rather than in winter.Winters probably were drier during the LIA due to anticyclonic weather patterns over the North Atlantic.
The higher number and thickness of sand layers in the Flugarth record in recent times confirmed that storm tracks shifted further north, as shown by studies summarised in the Intergovernmental Panel on Climate Change (Seneviratne et al., 2021).More severe storms and thus higher storm flooding can be expected in the North Atlantic realm for the future as the combination of rising ocean temperatures and cold Arctic air from the north could form intense low-pressure systems (e.g.Goslin and Clemmensen, 2017).

Figure 1 .
Figure 1.Map of the study area.(a) Three different storm track types cause storm surges along the North Sea coast: the Scandinavian type (yellow), the Skagerrak type (blue) and the Jutland type (pink) (from Wöffler, 2016, modified after Petersen and Rohde, 1991).(b) The Shetland Islands with the study area in north-western Mainland.(c) Sand Voe bay with yellow areas showing the calcareous fixed dune mapped by Dargie (1998).(d) Location of Loch Flugarth and position of the sediment cores (FLUG 1-4), and the modern environmental samples M1 and M2.[Color figure can be viewed at wileyonlinelibrary.com]

Figure 2 .
Figure 2. (a) rBacon age-depth model.All considered radiocarbon dates are in blue, while the non-radiocarbon ( 137 Cs) age-depth information is in green.The posterior age-depth model with the median (red dashed line) and the 2σ range runs from bottom left to top right.The horizontal bars in the figure represent event deposits (slumps).(b) rBacon age-depth model of the upper 15 cm with the non-radiocarbon ( 137 Cs) information only.(c) Vertical 137 Cs activity distribution in sediment core FLUG 2. Measured activity concentration in bright blue columns; modelled activity concentration in black diamonds.Time markers are presented by green horizontal lines.[Color figure can be viewed at wileyonlinelibrary.com]

Figure 4 .
Figure 4. Sedimentological results of FLUG 3. From left to right: photography, CT scan, lithostratigraphy (Log), density (blue line), CIE colour *b covering the blue-yellow axis colour spectra (yellow line), fractional porosity (pink line), and a heatmap showing the grain-size distribution frequency and the mean grain size (black line).[Color figure can be viewed at wileyonlinelibrary.com]

Figure 5 .
Figure 5. End-member modelling analysis showing: (a) the loadings of grain-size distribution frequency with the three different end-members (EM) together with corresponding samples and (b) the down-core scores for each end-member.[Color figure can be viewed at wileyonlinelibrary.com]

Figure 6 .
Figure 6.Factor map showing the results of the hierarchical clustering of principal components (HCPC).(a) Factor map of the HCPC with four different clusters.Numbers indicate the depths (individuals) in cm b.s.(b) Core FLUG 3 with provenance of samples shown in the factor map. Colours of the frames refer to colours and symbols of grouping in the HCPC plot.[Color figure can be viewed at wileyonlinelibrary.com]

Figure 8 .Figure 9 .
Figure 8. Overview of the periods with high storm frequency (red bars), age of event deposits (black points indicate the mean age), the main storm events, and the number of storms per decade and per half-century after Lamb (1991) between 1500 and 1989.[Color figure can be viewed at wileyonlinelibrary.com]

Figure 10 .
Figure 10.Comparison between NAO, different storm studies and the data of the current study over the last ca.1500 years (modified after Orme et al., 2016).Larger climate anomalies such as the MWP and the LIA are presented in red and blue, respectively.(a) Sedimentological signal of storm flooding by EM1 (current study).(b) Geochemical signal represented by XRF Si/Ti ratio (current study).(c) Number of storms per 20 years from historical documents for northern Scotland and the northern North Sea(Lamb, 1991).(d) Number of studies reporting sand transport in Scotland and Northern Ireland in the respective year(Orme et al., 2016).(e) Cliff-top storm deposits on the Shetland Islands(Hansom and Hall, 2009).(f) Winter (DJF) NAO reconstruction(Trouet et al., 2009).(g) Smoothed Greenland sea salt influences (ssNa+) of core GISP 2(Meeker and Mayewski, 2002).(h) Stalagmite proxy from a Scottish cave(Baker et al., 2015).(i) Reconstruction of Atlantic water inflow in the Norwegian Sea derived from concentrations of the coccolith Gephyrocapsa muellerae(Giraudeau et al., 2010).[Color figure can be viewed at wileyonlinelibrary.com]

Table 1 .
14C data of samples from the current study.
*represents a "negative" post-bomb age.†outlier,not used for the Bayesian age-depth model.
, and thus indicate organic-rich, peaty layers in the current study.Since Ca is often associated with carbonates and is used as an indicator of marine origin, it could either indicate © 2023 The Authors Journal of Quaternary Science Published by John Wiley & Sons Ltd.J. Quaternary Sci., 1-17 (2023) [Color figure can be viewed at wileyonlinelibrary.com] 1500 YEARS OF NORTH ATLANTIC STORM FLOODING 11 10991417, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/jqs.3568by NHS Education for Scotland NES, Edinburgh Central Office, Wiley Online Library on [24/09/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License © 2023 The Authors Journal of Quaternary Science Published by John Wiley & Sons Ltd.