Resolving tsunami wave dynamics: Integrating sedimentology and numerical modelling

Tsunamis are a major hazard along many of the world's coastlines. To understand the impact of these events, a sufficiently long record of previous events is needed, which can be provided by their sedimentary deposits. A number of past events have left extensive sedimentary deposits that can be used to understand the hydrodynamics of the tsunami. The ca 8.15 ka Storegga submarine slide was a large, tsunamigenic mass movement off the coast of Norway. The resulting tsunami had estimated run‐up heights of around 10 to 20 m on the Norwegian coast, over 30 m in Shetland and 3 to 6 m on the Scottish mainland coast. New cores were taken from the Ythan Valley in North‐East Scotland, where Storegga tsunami deposits have previously been found. High‐resolution sedimentary analyses of the cores, combined with statistical (changepoint) analysis, shows signatures of multiple waves. Moreover, detailed CT scans of the erosional basal surface reveal sole marks called skim marks. Taken in conjunction with the grain size and sedimentary fabric characteristics of the tsunami deposits, this indicates that the flow exhibited a high‐concentration basal component, with an initial semi‐cohesive phase and that deposition was dominantly capacity driven. A multiple wave hypothesis is tested by creating a high‐resolution numerical model (metre‐scale) of the wave inundation, coupled to a previously published regional model. The inundation model confirms that multiple waves passed over the site in agreement with the sedimentological analysis. The sensitivity of the model to the reconstructed palaeocoastal geomorphology is quantitatively explored. It is concluded that local palaeogeomorphological reconstruction is key to understanding the hydrodynamics of a tsunami wave group in relation to its sedimentary deposit. Combining sedimentological data with high‐resolution inundation modelling is a powerful tool to help interpret the sedimentary record of tsunami events and hence to improve knowledge of their risks.


| INTRODUCTION
Tsunamis are rare events (Scheffers & Kelletat, 2003).In order to understand these events, sedimentary deposits from palaeotsunamis can be studied to increase the time period over which extreme events can be assessed.The majority of studies on tsunami sedimentary deposits are either restricted to understanding the inundation distance or maximum run-up height (Long et al., 2016) or to interpreting the depositional characteristics of a tsunami deposit in terms of flow dynamics (Jaffe et al., 2012).However, to fully understand how the sediment was deposited, and to gain a deeper insight into the hazard posed, the hydrodynamics of the wave(s) that deposited those sediments must also be understood (Naruse & Abe, 2017).This remains a challenge as it is not yet fully understand how sediments are deposited under tsunami waves (Tappin, 2007) as they are rare events and hence difficult to study.Multiple depositional mechanisms have been proposed but most observed sedimentary features could also be the result of deposition of storm events and therefore careful consideration of the depositional setting and application of fundamental process sedimentology is required (Shanmugam, 2012).
Numerical models are useful when studying the sediment deposition that occurs as a result of a tsunami.Those models can be forward models (Bosnic et al., 2021;Sugawara et al., 2014) or inverse models (Jaffe et al., 2012;Mitra et al., 2020;Spiske et al., 2010;Tehranirad et al., 2021) whereby flow quantities are derived from the sediments.The methods in the inverse models rely on either a necessarily simplified model of sediment transport and settling (Jaffe et al., 2012;Tehranirad et al., 2021) or machine learning to obtain data from a large database of hypothetical deposits (Mitra et al., 2020).However, linking both forward and inverse models can give greater insights into the resulting sedimentary deposition than a single approach only (Bosnic et al., 2021).Regardless of the type of approach taken, a key part of any interpretation is a thorough understanding of the depositional mechanisms that occur during tsunami inundation.To date, this is poorly known (Costa & Andrade, 2020), which remains a challenge when interpreting palaeotsunami deposits (Yamamoto et al., 2021).
The Holocene Storegga tsunami, which occurred around 8150 years ago (Bondevik et al., 2003), was generated by a large submarine slide (Bondevik et al., 1997).The wave spread across the Norwegian-Greenland sea (Bondevik et al., 2005;Haflidason et al., 2005;Harbitz, 1992), affecting a region of 95,000 km 2 (Bondevik et al., 2005).Many tsunami deposits from the Storegga slide-generated wave have been found across the region, including Scotland (Dawson et al., 1988;Dawson & Smith, 2000;Long et al., 1989Long et al., , 2016;;Smith et al., 2004;Tooley & Smith, 2005), northern England (Boomer et al., 2007), Norway (Bondevik et al., 2003;Svendsen & Mangerud, 1990;Vasskog et al., 2013), Faroe Islands (Grauert et al., 2001) and Greenland (Wagner et al., 2007).Run-up heights are estimated to be over 30 m in some locations, particularly where the tsunami wave propagated large distances along Norwegian fjords (Vasskog et al., 2013).Run-up heights are estimated from sedimentary deposits along these coastlines by measuring the altitude of the deposit, which is then corrected for relative sea level (RSL) change over the last 8 kyr (Dawson et al., 2020).The deposits have been found using a variety of methods, including via lake coring (Rasmussen et al., 2018), sections (Dawson et al., 2006) and terrestrial cores (Smith et al., 1999).To date, no Storegga deposit has been analysed with the aim of establishing the hydrodynamics of the wave that generated the deposit by combining the sedimentology with numerical models.Moreover, previous attempts to understand the hydrodynamics via modelling of the Storegga tsunami waves have focussed on regional-scale models (Harbitz, 1992;Hill et al., 2014).
The first numerical simulations of the Storegga slide and the subsequent tsunami have shown how the wavetrain propagated, and have given reasonable matches to estimated run-up data, despite using present-day bathymetry (Bondevik et al., 2005;Harbitz, 1992).Later work extended these regional models using a multiresolution mesh and palaeobathymetry (i.e.present-day bathymetry corrected for RSL change), which improved the match to estimated run-up heights from sediments in some regions (Hill et al., 2014).All of these large-scale model runs did not include simulating the inundation of the wave onto the land (also known as wetting-anddrying) and all run-up heights were estimated from offshore wave height.A limited number of inundation studies have been attempted.For example Bondevik et al. (2005) coupled their regional model to a higher resolution inundation model to assess how the wave inundated an area of southern Norway.However, only one other inundation simulation to date has accounted for changes in coastal geomorphology and topography modelling, sedimentary analysis, Storegga, tsunami (Bateman et al., 2021), and hence most studies are limited in value for comparing to the sedimentary deposits found, and for understanding the hydrodynamics and the formation of the sedimentary deposit.Changes in coastal geomorphology may affect bathymetry, tidal conditions and coastal wave height, thus potentially having significant impacts on the accuracy of numerical models.These limitations, along with uncertainties in relative sea-level changes, may be a potential reason for why models of the Storegga tsunami do not fully agree with wave height estimates from sediment deposits (Dawson et al., 2020).
Since the Last Glacial Maximum there have been large changes in global mean sea level, which have been modified at the local and regional scale by changes in land height as the ice-sheets retreated (glacial isostatic adjustment, GIA; Shennan et al., 2018).These changes in RSL, together with sediment supply, control changes in coastal geomorphology.To fully reconstruct the inundation of a palaeotsunami, the land surface over which it flowed must also be reconstructed.Previous work on the Storegga tsunami has only accounted for the RSL change.Hill et al. (2017) simulated inundation across three regions of Scotland with a palaeobathymetry reconstructed from modern bathymetry combined with RSL change.They found that the model predicted very few of the coastal sites where the Storegga deposit was found were inundated unless sea level was higher at the time of the event than as assumed by their palaeobathymetric reconstruction.The level of relative sea-level at sites is a major driver of discrepancies between models and sedimentary data (Dawson et al., 2020).However, changes in coastal geomorphology and sediment accumulation were not accounted for, which presents a challenge to numerical modelling at sites where sedimentary deposits are found (e.g.dunes, salt marshes, tidal flats, etc.).In recent work Bateman et al. (2021) showed that correcting for geomorphological changes was key to attaining a reasonable match between a relatively coarse resolution inundation model and sedimentary data.
Here, new sedimentary cores are described from North-East Scotland, UK, that contained the Holocene Storegga tsunami deposit.Along with detailed analysis of the sediment, changepoint statistical analysis of grain size, is used to infer deposition mechanisms by three or four waves at each site from those sediment cores.High-resolution numerical modelling is then used, based on a variety of palaeobathymetric reconstructions, to understand how a wave group from the Storegga tsunami might have inundated the site.The results from numerical modelling are then integrated with information obtained from the sedimentary deposits to infer the depositional mechanism of tsunami deposits.

| SITE DETAILS
The Ythan estuary in North-East Scotland lies at the mouth of the Ythan River, where it drains into the North Sea (Figure 1).It is a mesotidal environment (mean tidal range of 3.1 m) and currently the tidal reach is 11 km inland (Stapleton & Pethick, 1996).A key feature of the modern estuary are the Sands of Forvie, which are a coastal sand dune system (Ritchie, 2000).The bedrock of the area comprises Dalradian metamorphic schist, overlain by glacial sands and gravels (Trewin et al., 1987).The modern estuary predominantly comprises salt marshes, tidal flats and coastal dunes (Sands of Forvie).The intertidal area was originally studied by Jamieson (1865) and then by Smith et al. (1983Smith et al. ( , 1999)), and summarised in Shennan et al. (2018) and Smith et al. (2018).These studies documented a RSL rise of ca 10 m from 10 ka, reaching a regional sea-level highstand ca 1.5 m above present at ca 7 to 5 ka.It is during the sea-level highstand that the Sands of Forvie are thought to have first formed (Ritchie, 2000).Lithostratigraphic surveys of the site reveal a succession of glacial deposits, overlain by peat horizons with transgressive surfaces of saltmarsh clays and tidal-flat silts (Smith et al., 1983(Smith et al., , 1999)).These deposits are widely overlain by a 5 to 7 cm thick sand horizon at elevations between 0.45 m and −0.2 m OD.Radiocarbon ages of peats are 8384 to 8019 cal yr BP (Smith et al., 1999) below the sand layer and 8175 to 7750 (Smith et al., 1999) cal yr BP above.Based on the ages of the well dated Green Moss deposits in Norway (Bondevik et al., 2012), these dates establish that the sand horizon was deposited by the Holocene Storegga tsunami.The Storegga deposit is described as a grey, micaceous, silty fine sand (Smith et al., 1983(Smith et al., , 1999)).At Waterside, the deposit forms a landward tapering wedge, which was previously interpreted as a storm surge deposit (Smith et al., 1983), but has since been re-evaluated as a tsunami deposit (Smith et al., 1999).Along Tarty Burn, the deposit is around 7 cm thick and has a 'normally graded grain-size profile', with a mean grain size of 70 μm at the base to 30 μm at the top of the deposit (Smith et al., 1999).

| Sedimentological analysis
Sections of the three cores extracted from the present day salt marsh using a Russian corer were used for analysis.Sections of each core (A004, A007 and B003) that contained deposits belonging to the Storegga tsunami event were initially visually identified.Core A007 was geographically oriented by marking north on the corer prior to rotation and extraction.The sections of each core containing the Storegga deposit were checked for any sedimentary structures or features.Images of the cores were recorded using a Panasonic LUMIX DMC-TZ60 digital camera.Subsequently, core sections were examined using a Zeiss Axio Zoom V16 microscope with an Axiocam 105 colour camera attached.
The density distribution of sediments through core A007 was investigated using X-ray computed tomography (CT).A Brivo CT385 medical CT scanner (GE Healthcare) was operated at an X-ray voltage of 120 kV, tube current of 80 mA and focal spot size of 1.2 mm for an exposure time of 1 s.This enabled the reconstruction of 96 mm diameter, 110 mm deep sub-samples of the core with a pixel resolution of 0.188 mm and a slice separation (or voxel depth) of 0.625 mm.This resolution is high enough to see sedimentary features in detail, but not sufficient to see individual grains.Image processing was conducted using Matlab v2020a.
Sampling of cores A004 and B003 began at ca 30 mm beneath the observable base of the deposit and ceased approximately 30 mm above the observable upper contact of the deposit.Samples were taken every 2.5 mm.The total core sections analysed were around 100 mm in length.Once removed, core samples were isolated in glassware and exposed to hydrogen peroxide (H 2 O 2 ) in order to remove any organic matter via decomposition.Once decomposition had subsided, as indicated by a lack of an effervescence reaction, samples were transferred into centrifuge tubes and washed using deionised water in a Fisherbrand vortex mixer.In order to prevent coagulation of particles within the clay fraction, three drops of sodium hexametaphosphate ((NaPO 3 ) 6 ) were added to each centrifuge tube.Sections were then analysed for the 0.5-2000 μm fractions in a Malvern Mastersizer with a Hydro 2000MU dispersion unit.Grain-size statistics were calculated using the G2Sd R package (Fournier et al., 2014).
In order to assess if changes in grain size were statistically significant, changepoint analysis was employed using the R (R Core Team, 2020) package changepoint (Killick & Eckley, 2014).Changepoint analysis estimates the location at which the statistical properties of a sequence of observations change and the package used does this in a maximum likelihood framework.The framework minimises a cost function, C:   where the data are y, which contains an unknown number, m, of changepoints, .Here, f(m) is a penalty function to guard against overfitting.Measurements in both the variance and mean of the mean grain size vertically in the cores were used to identify significant changes in the mean grain size through the cores.The analysis used PELT algorithm and the MBIC penalty function.See Killick and Eckley (2014) for details.

| Numerical modelling
Numerical models that simulate coastal ocean regions can be used to predict the altered hydrodynamics caused by the tsunami wave.In this case Thetis (Kärnä et al., 2018), a (2D and 3D) flow solver for simulating such coastal flows (Angeloudis et al., 2018); Baker et al., 2020), implemented using the Firedrake finite element partial differential equation solver framework (Rathgeber et al., 2017) was used to simulate the tsunami.Thetis was configured to solve the non-conservative form of the non-linear shallow water equations: where is the water elevation, H d the total water depth, g the gravitational acceleration (9.81 m/s), the kinematic viscosity of the fluid and u the depth-averaged velocity vector.The Coriolis term is represented as f u ⊥ , where u ⊥ the velocity vector rotated counter-clockwise over 90 o .In turn, f = 2Ωsin with Ω corresponding to the angular frequency of the Earth's rotation and the latitude.Bed shear stress ( b ) effects are represented through the Manning's n formulation as: where n is the Manning's friction coefficient.The model is implemented using a discontinuous Galerkin finite element discretisation (DG-FEM), using the P 1DG − P 1DG velocity-elevation finite element pair.A fully implicit second-order DIRK22 time discretisation was used with a fixed time step Δt = 3 s to capture small temporal-scale wave dynamics on a mesh with minimum resolution of 20 m.Wetting-and-drying inundation processes are represented according to the formulation of Kärnä et al. (2011).Finally, the discretised equations are solved using a Newton non-linear solver algorithm using the PETSc library (Balay et al., 2017).Note the model does not include algorithms to simulate the bedload transport of sediment as this is challenging and requires calibration even when comparing to laboratory data (Yamamoto et al., 2021).
To simulate the Storegga tsunami inundation in the Ythan Valley a multi-scale mesh was constructed to encompass the valley and the coastline adjacent (Figure 2).The mesh was constructed using the 30 m topographic contour to form the landward boundary, which is then joined to an arc form on the seaward boundary.Boundaries and mesh creation were carried out in QGIS (QGIS.org,2020), using qmesh (Avdis et al., 2018) and Gmsh (Geuzaine & Remacle, 2009).Mesh resolution is a linear function of distance from boundary and water depth, with increased resolution in the region of interest.Resolution varied from 20 m in the area of interest around the Ythan Valley to 1 km away from the coastline in deeper water.Using the HRDS Python package (Hill, 2019), the full range of topographic and bathymetric data (see below for full description of these data) were dynamical loaded into Thetis flow solver during the simulation so that to ensure the finest resolution data is selected at each point in the mesh.The seaward boundary is forced using water elevation data from the regional Storegga simulation of Hill et al. (2014) via oneway coupling.The data from Hill et al. (2014) were interpolated onto the node points of the mesh via linear 2D interpolation in space.The wave height data were then linearly interpolated in time to the 3 s timestep of the Thetis flow solver.This is the same methodology as Bateman et al. (2021).
The initial bathymetric/topographic data are a mix of spatial resolutions.The topography data were derived from the 50 m Ordnance Survey DEMs and the Scottish Government's LIDAR data at 1 m resolution.Bathymetric data were derived from OceanWise data, which have a ca 30 m resolution.These data were adjusted to the same datum (Ordnance Datum).To adjust these data further to account for RSL and geomorphological changes five different scenarios were created to represent different levels of recreation of the palaeolandscape (Table 1).The same strategy as Hill et al. (2014) was employed to ensure consistency in the boundary conditions and the simulated scenarios presented here, which was to apply the GIA reconstruction of Bradley et al. (2011) to recreate the palaeobathymetry/ (1) topography.This simulation is labelled as the palaeo set-up (Table 1).A set of further palaeobathymetric-topographic reconstructions were created to account for the uncertainty in geomorphological changes not included in the palaeo reconstruction.First, the river aggradation that has occurred over the last 8150 years was removed by assuming the river does not move in a lateral direction, but the river mouth would be at sea level at 8.1 ka.A linear function of distance was used in the upstream direction of the coastline which varied from 1 at the 0 m contour to zero at 15 km away from the coast (which is outside the domain of interest):  where D is the distance from the 0 m contour (m), H the original palaeotopographic reconstruction (modern day topography minus RSL) in metres and m is: where R h is the height difference of the palaeoriver mouth to palaeosea level (m) and d is distance away from the centre of the river channel (m).This function simulates a river channel with the maximum removal of aggradation at the centre of the channel to take the channel down to sea level and reducing to zero at 500 m away from the centre of the channel.The resulting topographic/bathymetry is labelled as the palaeo-river set-up.The second modification was to remove the Sands of Forvie (SoF), which formed ca 6 ka (Ritchie, 2000).An approximate polygon was drawn over the modern day area and used to create a linear function from 1.0 inside the polygon to zero 500 m away from the outer edge of the polygon.This was multiplied by the height above sea level of the palaeo-river set-up within that polygon.The spit was removed in stages by subtracting 0.5, 0.8 and 1.0 of the total within-polygon height.The full reconstruction where the whole spit height is removed is labelled as palaeofull (Figure 3).Simulations where 0.5 (50%) or 0.8 (80%) of the spit removed are labelled as palaeo-0.5 or palaeo-0.8respectively.All simulations are summarised in Table 1.

| Sedimentology
Visual, CT and microscopic analysis conducted on the cores obtained from the Ythan estuary allowed the identification of several features and structures.Visual inspection of each Storegga deposit showed an overall normally graded grain-size profile from base to top.A recurrent feature throughout each tsunami core is the presence of vegetative and disarticulated shell material (Figure 4).Within the tsunami deposit in B003, the basal contact between the tsunami deposits and the underlying units appears sharp, with minor irregularities in the form of the contact (Figure 4).The basal contact to tsunami deposit A004 is also sharp, with marked erosional features.Within both cores, a change in compositional properties is reflected at this contact, with the underlying units richer in peat and devoid of siliciclastic material.The upper contact between the tsunami deposit and the overlying units whilst distinct, is not as sharp as the basal contact and is termed gradational.This contact, like the lower boundary, marks a compositional change between the two units, with the tsunami deposit being richer in siliciclastic material than the overlying unit.The tsunami deposits in cores A004 (Figure 5) and B003 (Figure 6) are characterised by silty sands and sandy silts, respectively.The tsunami deposit was determined to be 74.5 mm and 69.5 mm thick in cores A004 and B003, respectively.Both cores show similar grain-size changes.The base of the deposits are inverse graded over the lowermost 1-3 mm.Core A004 shows a further coarsening to 80% sand approximately 10 mm from the base.The grain size in core A004 then remains almost constant in terms of fine sand, but there is a steady decrease in medium sand from 58 mm below the top of the deposit to 24 mm below the top.The grain size then stays relatively constant until the top (Figure 5).In contrast, B003 shows a normal grading, with medium sand absent from 26 mm below the deposit top (Figure 6).However, there is then a clear inverse grading at 11 mm below the deposit top, followed by an abrupt normal grading to the top of the deposit.Neither core showed any clear sedimentary structures other than the erosional features at the base.
The cores show distinct changes in grain-size statistics below, within and above the tsunami deposit.Calculating sorting ( g ), skew (Sk g ) and kurtosis (K g ) using the geometric Method of Moments (Blott & Pye, 2001) shows a distinct shift from very poorly sorted below the deposit to poorly sorted deposition within the deposit (Figures 5  and 6).The kurtosis shifts from platy/meso-kurtic before and after the tsunami deposit to very leptokurtic within.The skew also shifts from largely symmetric to very fine skewed within the deposit.Similar changes are observed within both cores.
To analyse the observed grain-size changes further, the changepoint analysis shows each core had at least three significant changes from the base of the tsunami deposit to the identifiable top.Both cores show a change during the abrupt coarsening of grain size from the peat to the sand as expected.Additionally, both cores show three distinct levels in the mean/variance.In core B003, the first two of these are clear, but the final changepoint is driven by an abrupt and short-lived increase in mean grain size at 12 mm depth, rather than a visibly clear shift in mean grain size (Figure 6).There is also a very small section at 68.5 mm depth which may be of importance.Referring to the grain-size data there are, however, distinct increases in the proportion of sand where the changepoint analysis has indicated a statistically significant change.In contrast, Core A004 shows three clear signals, including the initial peak in grain size (Figure 5).Each statistically significant region has a clear flattening of mean grain size, followed by a decrease in grain size where the next changepoint is detected.

| CT scans
To determine the initial mode of sediment transport, the CT scan was focussed on the basal surface of the deposit in Core A007.The density variation in the CT scan can be plotted as colour range, with higher radiodensity representing sand (yellow) and lower radiodensity (green to blue) representing organic material (peat) or clay.The CT scan of the bottom 7 mm of the A007 core shows a clear increase in sand (Figure 7).There are several patches of sand at the northern edge of the deposit, which are elongate in a roughly west to east direction.The extreme high density on the western edge, below the basal surface, is probably an artefact from coring as it sits below the erosional base.Unlike the grain-size analysis of A004 and B003, the CT scan of A007 also show inclusions of peat clasts at heights of 9 mm, for example Figure 7A dotted circles, and more peat-rich layers at heights of 20 mm and around 45 mm indicated by the lower density values (Figure 7B).

| Numerical modelling
The palaeo-full numerical modelling shows the initial smaller positive wave producing a small increase in water level along the coastline, but this wave does not inundate the cores sites (see Video S5).This is the leading wave from the Storegga collapse, which is mostly directed towards Greenland (Hill et al., 2014).There is then a large drawdown of before the first tsunami wave then travels south along the coast, inundating the site.This is also the primary Storegga-induced wave, which according to Hill et al. ( 2014) is derived from the wave directed initially towards Norway.Subsequent waves, which are generated by reflection from the neighbouring coastlines and internally within the valley, then inundate the site.These waves are smaller than the first tsunami wave, but they result in a larger overall inundation depth as they impact already inundated areas.After 5 h of simulated time, the water has yet to fully recede.The flow patterns associated with the waves are relatively straightforward, with the main flow along the axis of the valley.There are small divergences of flow produced around and alongside of topographic highs (Figure 8).In detail these could affect the direction of flow of the wave, depending on the exact palaeotopographic reconstruction used, but the direction of the flow as the leading edge of the wave passes over the site is towards the north-west and turns west as the wave recedes.The results of the numerical models show a wide contrast in wave behaviour, depending on the palaeobathymetric reconstruction used.Using a simplistic palaeogeomorphological reconstruction (palaeo) results in no inundation of either core location, nor of the previous two studies where Storegga deposits are found (Video S1).Accounting for river aggradation (palaeo-river) does result in the site being inundated (Video S2).However, only two waves are observed that have sufficient velocity to move the grain sizes observed in the core (Figure 9).Once the removed, even partially, all three simulations showed three waves inundating the site (Videos S3-S5).Full removal of the spit and 80% of the spit removed (palaeo-full and palaeo-80) both give very similar results in terms of the inundation depth, velocities and number of waves present (Figure 9).Removal of the SoF is reasonable when considering similar coastal features along the eastern Scottish coast appeared to have form during the mid-late Holocene, post-dating the Storegga tsunami (Firth et al., 1995).Three clear waves are recorded at the site, with an additional fourth wave indicated by a peak in velocity observed at around 12,500 s.Each inundation is larger than the previous wave as the subsequent waves occur before water has fully receded.Flow speeds at the site for these two simulations show an initial peak in velocity associated with the incoming wave, followed by smaller subsequent peaks for the subsequent inundation events.However, the narrow section of the valley to the south of the site, accelerates flow to over 2.5 m/s during the third wave (Figure 8); around three times the value at the core site.The palaeo-0.5 simulation also shows three inundation periods at the site, but these are less well defined, particularly in the flow depth, with no clear reduction in depth between waves to allow for deposition of sediment F I G U R E 9 Flow depth (A), flow speed (B), and bed shear stress (C) at the core site.C has grain sizes shown at the critical bed shear stress (VFS, very fine sand; MS, medium sand; CS, coarse sand; VCS, very coarse sand; VFG, very fine gravel).Depending on spit height the maximum flow depth varies from 2.5 m (where the spit is fully removed) and just under 1 m (where the SoF is fully intact).Peak flow speeds are 1.2 to 0.4 m/s.The speed generally falls with increasing wave height, meaning the first wave, though not the largest in terms of flow depth, usually shows the fastest flow speed.This translates to the flow being capable of transporting very fine gravel for the highest velocity simulated.from suspension.between the simulations and the effect the spit can be seen in Video S6.
The bed shear stress as defined by Equation ( 4) can be converted into maximum grain size that can be mobilised and eroded from the source region.Bed shear stress is sufficiently high to move all grain sizes found in the two sites under either suspended load or bedload conditions, or a combination of both.It is therefore clear that the grains found are not transport dependent, but source dependent.Each successive wave has a lower bed shear stress as a result of both decreasing velocity and increasing water depth.The direction is such that the flow is to the north, but then turns westwards just south of the deposit site, and towards the east, north of the site (see Video S7 and Figure 8).

| DISCUSSION
High-resolution analysis of sedimentary deposits dated to the time of the Storegga tsunami indicate that multiple waves inundated the site.There is clear evidence of at least three waves via changepoint analysis at both sites, with some evidence of up to four waves in core B003.The CT scans, sedimentological analysis and numerical model are broadly consistent with a clear depositional mechanism for the sediments.

| Interpretation of the sedimentary deposit
The two cores where grain size was analysed (A004 and B003) show an erosional base, and in addition, orientated linear erosional marks were identified on the CT scan of core A007.These are interpreted as tool marks, specifically skim marks, formed by larger grains that would occasionally arc downwards and erode the bed.Skim marks are associated with high-concentration semi-cohesive flows that are able to provide buoyant support to the clasts (Baas et al., 2021;Peakall et al., 2020).This erosional surface and associated tool marks are probably formed by the front of the flow which will preferentially erode and concentrate the most mobile sediment (Baas et al., 2021).Observations from videos of the 2011 Tohoku-oki tsunami, Japan, traversing flat fields, show the front of the initial wave transforming into a high-concentration debritic head as a result of large-scale erosion and entrainment of fines (Tappin et al., 2012).Such flow transformation through incorporation of fines is well-known in subaqueous sediment gravity flows (Talling et al., 2004), producing either fully laminar flows or flows transitional in their behaviour between laminar and turbulent states (Baas et al., 2009) that are linked to the formation of these tool marks (Baas et al., 2021;Peakall et al., 2020).The deposits of such a high-concentration head, whilst not observed here, would be expected to be erosionally based, increasingly poorly sorted upwards, clast-rich and to pinchout abruptly (cf.Bondevik et al., 2003;Costa et al., 2021, their Figure 5G,H,I).
The base of the deposits overlying the erosional surface are inverse graded over heights of 1 to 3 mm, representing either an increase in flow velocity (Mulder et al., 2003) or deposition from laminar shear layers in decelerating highconcentration basal layers (Sumner et al., 2008).The overlying deposits show some differences in grain-size profiles between the two cores, despite being collected less than 200 m apart.Both cores show poor to very poor sorting, with grain size varying from silt to medium sand and a general decrease in mean grain size towards the top of the deposit.The very leptokurtic distribution shows that the distribution is far more peaky than a normal distribution, reflecting that the centre of the distribution is much better sorted than would be expected in a normal distribution, indicating a mixing of sources (Folk & Ward, 1957).This is interpreted as evidence of erosion and entrainment of beach deposit or very well sorted shallow marine sand, which then incorporated estuarine/river sediment rich in fine-grained material as the wave inundated the Ythan Valley.If the Folk and Ward (1957) graphical method is used then these sediments are estimated as well sorted, which is misleading as the combination of silts and sands suggests that sorting is poor, as confirmed by the Method of Moments.This graphical method under-represents the tails of the distribution and may be particularly misleading in tsunami sediments such as this, where kurtosis is very high (Swan et al., 1978).Core A004 shows clear evidence of layers within the grain-size data, although not visually.There are three layers in total, and each layer corresponds to a decrease in mean grain size of around 20 μm in each layer.In contrast, B003 shows a more gradual decrease with only a short section of core that has a clear constant mean grain size.The presence of a wide range of grain sizes from silts to medium sands, combined with the sorting and kurtosis of the deposit, suggests that deposition took place through a loss of capacity of the flow to carry the total sediment load, thus deposition occurs at velocities where a clear current with the same velocity would be able to suspend the sediments (Hiscott, 1994).The basal layers, and in the case of core A004 the upper layers, show a near-constant grain-size distribution, suggesting that deposition took place from a high-concentration (>30-45 vol.%) flow characterised by hindered settling (Cartigny et al., 2013;Druitt, 1995).This evidence for high-concentration suspensions is in accordance with the inverse grading at the base of the bed, and the absence of tractional structures.Higher up, B003 is characterised by progressive signatures, and thus some grain-sorting occurred, however, deposition still covers a wide range of grain sizes, thus suggesting that capacitydriven sedimentation of a suspension was still occurring.Such deposits are associated with lower, but nonetheless still high, concentrations.Modelled basal shear stress values, associated with the three tsunami waves, are based on clear water flows and predict that all grain sizes should be in suspension during most of the duration of the three waves.The velocity simulated at the core site is large enough to move very fine gravel and the backwash is fast enough to move only very fine sand, based on the criteria for bedload transport.There is therefore time for deposition of sediment between each inundation, but these simulations suggest, based as they are on clear-water flows, that settling of the very fine sand and silt grains would not be expected.The finer components would only settle after the three waves have passed.The sedimentological observations supporting primarily capacity-based deposition are therefore in agreement with these hydrodynamic outputs, demonstrating that deposition can take place under tsunami waves at appreciable velocities (tens of cm/s).Such capacity-driven sedimentation in tsunamis has also been proposed by Naruse and Abe (2017), albeit they argue that this occurs at lower concentrations.These sedimentation patterns also bear similarities with deposition from many suspension-driven turbidity currents where capacitydriven sedimentation leads to a wide range of grain sizes being deposited at velocities that a clear-water flow would maintain in suspension, with the stratification of different grain sizes leading to normally graded sequences as a flow decelerates (Hiscott, 1994).
The changepoint analysis is interpreted as representing the number of waves that inundated the site assuming preservation of deposition under all waves; either settling of sediment in stages or as the impact of three or four waves at the site in terms of pulses of (re)suspension and deposition.The reducing grain size in each 'pulse' can be interpreted as decreasing wave energy over time.The overall interpretation of the sedimentary succession in both cores is broken into two main phases.First, a high-concentration, semi-cohesive, flow front cut the initial erosional surface and associated tool marks.Second, deposition of source sediment and associated predominantly finer-grained erosion products, took place from a high-concentration basal layer at some point behind the front of the initial tsunami wave.It is envisaged the whole water column to be carrying sediment at this time, but to exhibit prominent stratification.This initial deposit would have been laid down on the incoming wave but the upper section was probably reworked as this wave started to recede.The next two waves (reflections) then inundated the site.These have a lower bed shear stress (due to both lower velocity and deeper water) and subsequent units are interpreted to be deposited either as sheets of sediment carried towards the bottom of the water column, or where grading is observed, settling from suspension.In core A004, the sheet-like deposition dominates, producing the clear steps in mean grain size observed (Figure 5).Core B003 shows a short section with near constant mean grain size, which can be interpreted as sheet-like deposition, followed by a gradual decrease in mean grain size which can be interpreted as deposition from suspension.Each subsequent wave deposited a slightly finer overall grain size which may indicate reworking of the initial deposit, rather than a fresh input of sediment from the source region.The visually sharp top to the deposit and the absence of a finer cap may indicate that the wave train passing the sites created a 'vertical pumping' effect (O'Hara Murray et al., 2012), which would keep finer-grained sediments in suspension, preventing them from settling.These would then be transported seaward from the site as the wave receded.This corresponds to observations from the 2011 Japan tsunami where the backwash following the tsunami carried a large amount of suspended sediment offshore (Udo et al., 2015).The interpretation presented here fits with the idea of the waves passing over the site being largely reflective in nature, rather than large waves as part of the tsunami wave train.There are subsequent waves from the main tsunami wave train, but at the point at which it arrives at the site the waters have receded enough to disconnect the valley from the ocean and the subsequent wave is not large enough to cause further inundation (end of Video S5).The differences in the two cores is probably due to very localised geomorphological differences at the time of deposition, despite the close proximity.That localised deposition could be a result of localised scouring by the initial wave.However, this is not recreated in the numerical model as both locations present almost identical inundation and velocity histories, but the model does not simulate scouring or other processes that may create such localised changes in bed topography.
Of particular note here are the sole structures imaged via CT scans.Grooves (continuous on the scale of the outcrop) have been observed under putative Pleistocene tsunami deposits (Madeira et al., 2020), and 'megatsunami' deposits associated with the KT mass extinction (Smit et al., 1992(Smit et al., , 1996)).Extensive planform bedding plane exposures at the base of the Cretaceous-Palaeogene megatsunami sites reveal a range of discontinuous tool marks including skim marks up to several metres in length, and possible prod marks (Olsson et al., 1996).The marks are referred to as grooves by Olsson et al. (1996) but the planform maps clearly show them to be discontinuous features that are largely symmetrical in longitudinal form, and thus are better termed skim marks 1982; Peakall et al., 2020).The lying tsunami beds are associated with rip-up clasts, where eroding cohesive substrates (Madeira et al., 2020;Olsson et al., 1996).Whilst bedding plane surfaces are not well exposed in modern tsunami deposits, rip-up clasts are frequently observed (Long et al., 2016;Peters & Jaffe, 2010).The data presented here form the first evidence for sedimentary sole structures found in Holocene tsunami deposits to the authors' knowledge.Those are interpreted from Figure 7A as bright high density linear structures (circled) with rip up-clasts indicated as low-density areas (dotted circles).As discussed earlier, these features are associated with high-concentration, transitional, semi-cohesive flows, during an initial erosional phase (Peakall et al., 2020), and thus most probably form under a debritic flow head.The orientation of these features is largely east-west in core A007 indicating flow in towards either the east or west.The numerical modelling (palaeo-full) indicates flow initially to the north-west, before turning west as the wave starts to recede.Whilst not in total agreement, minor changes to the topography over the last 8000 years could result in the numerical model producing flow in these directions on the initial wave.

| Impact of palaeogeomorphological reconstruction
In previous studies, two other sites in the Ythan Valley have been found that record the Storegga deposit (Smith et al., 1983(Smith et al., , 1999)).These are Waterside and Tarty Burn.These sites were not included when reconstructing the palaeobathymetry/topography, so they act as independent checks on the accuracy of this reconstruction.In the simulations the maximum elevation of the wave almost inundates these two sites.At Waterside, there is a steep side to the valley, which may be a result of the reconstruction process and cosine function used to de-aggrade the river.The site at Tarty Burn is inundated but only to a very shallow (centimetre) depth and bed shear stresses recorded there (not shown) indicate the water velocity was not high enough to move the grains found in the deposit.This only occurs in the model when the spit is removed entirely (see Video S7).It is clear from Figure 10 that the reconstruction still contains some modern elements of topography.The curve just north of Waterside is a bridge, which has not been removed fully during processing of LIDAR data to create a DEM.This will have some impact on the wave around the Waterside site, as F I G U R E 1 0 Three dimensional representation of the maximum wave height from model palaeo-full.The view is looking north-east.All data are exaggerated by a factor of 10 in the vertical.The maximum wave height shows a general decrease as it travels up the Ythan Valley.The sites described here are labelled, with the previous locations of where a Storegga deposit was found (Tarty Burn and Waterside) shown using a white cone.The site described here is shown with the blue cone.The location is denoted by the tip of the cone.
there is change in maximum wave height at this location.The impact is negligible for the core sites in this study.Clearly, further work is required to optimise the reconstruction such that the model can accurately represent the local-scale deposit found at all three sites in the Ythan Valley.

| Modelling uncertainties
The uncertainties in modelling the Holocene Storegga tsunami come from a number of inter-related sources.First, the source wave is complex and heavily dependent on the type of mass flow and how it fails (Løvholt et al., 2017).The source wave used in this study is from Hill et al. (2014) who used a relatively simple solid block model.Subsequent work has shown that modelling a deformable slide, even in the same numerical framework as the Hill et al. (2014) study, can produce a more complex wave form (Smith et al., 2016).How this then translates to inundation at the high resolutions used here is an open question.Second, the timing of the tsunami is important when attempting to match to sediment deposits, particularly in meso-tidal and macro-tidal regions.The modern tidal range in the Ythan Valley is around 3 m, so even assuming no change in tidal range over the last 8 kyr, the tsunami run-up height could increase by up to 1.5 m assuming a simple relationship.It is know that the Storegga mass flow and tsunami occurred in the autumn (Rydgren & Bondevik, 2015), but the timing of the slide with respect to the tidal cycle is not known.The interaction of these two waves can affect both run-up height (Weisz & Winter, 2005) and current speed (Kowalik et al., 2006), which in turn will affect bed shear stress.Additional complexities of the discharge in the river at the time of the tsunami arrival were also ignored in this study.Moreover, the Manning's drag coefficient plays an important role in calculating bed shear stress and the water velocities in the wave.Here, a constant drag coefficient suitable for clean sand has been used, but vegetation and different bed characteristics will alter that.More work is needed to quantify the impact tidal dynamics and drag coefficient variability would have on the run-up of the Storegga tsunami in a tidal estuary, such as the Ythan Valley.Finally, the model here is non-dispersive so will not capture phenomena such as bore waves breaking (Madsen et al., 2008) which may increase tsunami amplitude during inundation.These sources of uncertainty, coupled with that from the palaeotopographic reconstruction means that a numerical model can never be a true reconstruction.However, it has been shown how using a numerical model can help inform the depositional mode of the tsunami sediment and help interpret grain-size variations through a deposit.This is a powerful combination of techniques to use to understand tsunami-related sedimentation.

| CONCLUSIONS
High-resolution numerical modelling can be a powerful tool in helping to understand sediment transport and depositional processes under tsunami waves, including in complicated topographic/bathymetric configurations.Here, it has been shown that multiple waves can be detected from high-resolution grain-size analysis of tsunami sedimentary deposits and statistical analysis of those data can be used to interpret associated waves of sediment transport.Detailed analysis using CT scans and sediment cores indicate that transport of sediment during this tsunami occurred as a sheet of high-concentration sediment, that initially is semi-cohesive, followed by dominantly capacity-driven deposition from this high-concentration flow.Depending on local conditions this may not be preserved, and can lead to different interpretations from deposits that are close to each other.When coupled with the numerical modelling a picture can be constructed of the whole wave, including when sediment is probably to be deposited, which is predominantly on the incoming wave, rather than the backwash.The modelling and sedimentological data show an excellent agreement once palaeogeomorphological changes are accounted for.These interpretations can be applied to other palaeotsunami deposits.However, key to the modelling is the palaeogeomorphological reconstruction.Whilst the model is not overly sensitive to changes in terms of broad scale interpretations, fine-scale interpretations are only possible when there is confidence in the reconstruction, including uncertainties.

ACKNO WLE DGE MENTS
This project formed part of LH's MSc project supervised by GR, JH and RG.We acknowledge funding awarded to JH, JP, GK, NB, EB, DH and RG from the White Rose University Consortium Fund (INUNDATION project).This project was undertaken on the Viking Cluster, which is a high-performance compute facility provided by the University of York.We are grateful for computational support from the University of York High Performance Computing service, Viking and the Research Computing team.We thank Prof. David Tappin for a constructive review.

FUNDING INFORMATION
White Rose Collaboration Grant: INUNDATION.

F
Location of the Ythan Valley, with satellite image (ESRI) showing the field locations in this work marked in purple.The locations of previous studies are indicated in blue (Tarty Burn and Waterside), along with the Sands of Forvie.The red box highlights the area covered below.Bottom right shows the locations of the three cores analysed in this work with environments labelled.

F
I G U R E 2 Computational modelling domain used in the study.(A) Black lines indicate a 'no flow' boundary, which was derived from the 20 m contour.The red box highlights the extent of map B. The red line indicates the domain boundary forced from the model of Hill et al. (2014).(B) Close up of the computational mesh using the Ythan Valley area with colours showing the fully reconstructed bathymetrytopography.The red box highlights the extent of map C. Mesh resolution is 20 m here with DEM resolution of 10-25 m.The location of the cores are labelled.(C) A close up of the field site where the high mesh resolution can be seen.Bathymetry c ○ British Crown and OceanWise, 2016.All rights reserved.Licence No. EK001-20180802.Not to be used for Navigation.Topographic data derived from Ordnance Survey data c○ Crown copyright 2020 and Scottish Government data.
Localised bathymetric and topographic reconstruction of the Ythan Valley.(A) The palaeo reconstruction, which is the modern bathymetry/topography corrected for relative sea level.(B) The full reconstruction after accounting for river aggradation and spit formation (palaeo-full).Both panels use the same spatial and colour scales.(C) The difference between (A) and (B).Most of the difference is located around the SoF, with smaller differences to account for river aggradation.F I G U R E 4 Samples of photographs taken of the two tsunami sediment cores.(A) Example of disarticulated shell from core A004.Scale bar is 500 μm.(B) Plant fragment from core B003.Scale bar is 1 mm.(C) Erosional base from core A004 highlighted by dotted line.There is clear erosion into the lower unit with a sharp contact.Scale bar is 500 μm.(D) Base of core B004.A sharp contact between the tsunami deposit and the underlying peat can be clearly seen, but lacks some of the erosional features observed in core A004.Scale bar is 500 μm.

F
Grain size and statistics of cores A004.Darker colours indicate coarser grain sizes (CS, coarse sand; MS, medium sand; FS, fine sand).The sedimentological statistics are sorting ( g ), skew (Sk g ) and kurtosis (K g ).Sorting: PS, poorly sorted; VPS, very poorly sorted.Skew: VFS, very fine skew; FS, fine skew; S, symmetric.Kurtosis: P, platykurtic; M, mesokurtic; L, leptokurtic; VL, very leptokurtic.The erosional base of the tsunami deposit is indicated with the red line and the upper surface by the yellow line.Depths are given with the top of the tsunami deposit as zero.

F
Grain size and statistics of cores B003.Darker colours indicate coarser grain sizes (CS, coarse sand; MS, medium sand; FS, fine sand).The sedimentological statistics are sorting ( g ), skew (Sk g ) and kurtosis (K g ).Sorting: PS, poorly sorted; VPS, very poorly sorted.Skew: VFS, very fine skew; FS, fine skew; S, symmetric.Kurtosis: P, platykurtic; M, mesokurtic; L, leptokurtic; VL, very leptokurtic.The erosional base of the tsunami deposit is indicated with the red line and the upper surface by the yellow line.Depths are given with the top of the tsunami deposit as zero.F I U R E 7 Horizontal slices and 3D orthoslices (B) of the CT data from core A007.The location of the slices in (A) are highlighted onto the axis in (B).The top of the deposit is very clear in (B).The base exhibits substantial erosional features.The circled areas in (A), towards the base, show examples of high density, linear features, with elongation in the east-west direction.Low density inclusions are also highlighted.HU refers to Hounsfield Units which are units of radiodensity, translated herein to grain size.F I G U R E 8 Instantaneous flow velocity vectors during the first (9660 s), second (10,740 s), third (12,060 s) and potential fourth (12,360 s) waves from model palaeo-full.The site of the three cores is shown by the magenta circle.Time in the top left corner shows time in seconds since simulation started.There is a clear divergence of flow in the vicinity of the site, the exact nature of which will depend on the palaeobathymetric/ topographic reconstruction.