Highly variable recurrence of tsunamis in the 7,400 years before the 2004 Indian Ocean tsunami

The devastating 2004 Indian Ocean tsunami caught millions of coastal residents and the scientific community off-guard. Subsequent research in the Indian Ocean basin has identified prehistoric tsunamis, but the timing and recurrence intervals of such events are uncertain. Here we present an extraordinary 7,400 year stratigraphic sequence of prehistoric tsunami deposits from a coastal cave in Aceh, Indonesia. This record demonstrates that at least 11 prehistoric tsunamis struck the Aceh coast between 7,400 and 2,900 years ago. The average time period between tsunamis is about 450 years with intervals ranging from a long, dormant period of over 2,000 years, to multiple tsunamis within the span of a century. Although there is evidence that the likelihood of another tsunamigenic earthquake in Aceh province is high, these variable recurrence intervals suggest that long dormant periods may follow Sunda megathrust ruptures as large as that of the 2004 Indian Ocean tsunami.

P rojections of fatalities due to catastrophic earthquakes and tsunamis will likely exceed 2 million lives in the twenty-first century 1 . Advances in geodesy and seismology have contributed to our understanding of rupture patterns of large earthquakes, but the devastation caused by the 2011 Tohoku-oki and the 2004 Indian Ocean tsunamis make it clear that estimates of earthquake size and tsunami potential are woefully inadequate. The repeat times of such giant tsunamis can occur centuries to millennia apart [2][3][4][5] and are not fully captured in historical and instrumental records 2,3 . A more refined understanding of the long-term variations in timing and recurrence of giant tsunamis is essential for producing realistic vulnerability assessments for coastal communities.
The great Sumatra-Andaman earthquake triggered a tsunami that devastated south and southeast Asia 5,6 . At the time, there was no known historic precedent for the 1,500 km rupture of the Sunda megathrust 5 , with slip exceeding over 20 m (refs 6,7). In the decade since the Indian Ocean tsunami, the search for prehistoric estimates of earthquake recurrence and tsunami potential remains elusive. Most reconstructions of past tsunami inundation are based on identifying anomalous beds of sand in low-energy environments, such as salt and freshwater marshes, coastal lakes or swales 8,9 . Prehistoric tsunamis have been identified using such geological records from northern Sumatra [10][11][12] , Thailand [13][14][15][16][17] , Andaman Islands 18 , Sri Lanka 19 , Eastern India 20 and the Maldives 21 , but the timeline of their reconstructions is limited or fragmentary, hindered by preservation problems, reworking and a lack of accommodation space 22 .
We identify coastal caves as a new depositional environment for reconstructing tsunami records and present a 5,000 year record of continuous tsunami deposits from a coastal cave in Sumatra, Indonesia (Fig. 1), which shows the irregular recurrence of 11 tsunamis between 7,400 and 2,900 years BP. The sedimentary record in the cave shows that ruptures of the Sunda megathrust vary between large (which generated the 2004 Indian Ocean tsunami) and smaller slip failures. The chronology of events suggests the recurrence of multiple smaller tsunamis within relatively short time periods, interrupted by long periods of strain accumulation followed by giant tsunamis. The data demonstrates that the 2004 tsunami was just the latest in a sequence of devastating tsunamis stretching back to at least the early Holocene and suggests a high likelihood for future tsunamis in the Indian Ocean. The sediments preserved in the costal cave provide a unique opportunity to refine our understanding of the behaviour of the Sunda megathrust, as well as study in detail the sedimentology and hydrological characteristics of tsunami deposits.

Results
Geologic setting. The coastal cave site is located along the northwestern coast of Aceh Province near the village of Lhong, 35 km south of Banda Aceh (Fig. 2). This segment of the Sunda megathrust ( Fig. 1) slipped as much as 20 m during the 2004 rupture 6,7 and produced nearly 1 m of subsidence. The 2004 tsunami inundated the cave and removed vegetation off the very steep limestone cliff to a height of B24 m above mean tidal level (MTL) which was over 10 m above the top of the cave entrance (Fig. 2). The cave entrance is 100 m back from the swash zone with a rock sill at its entrance that sits 1 m above mean tidal level (Fig. 2). The cave extends nearly 120 m into the cliff. We excavated six trenches at the rear of the cave (Fig. 2) and found sedimentary sequences up to 2 m thick above a limestone basement.  Tables 1 and 2). Basal rip-up clasts, lenticular laminations and fragments of weathered cave chalk are common in the 2004 sand bed in all trenches. The 2004 sand bed contains abundant, pristine foraminifera, mostly of benthic subtidal origin 23 , but with a notable planktonic offshore presence. Organic debris, transported into the cave by the tsunami and guano from the insect-feeding bats (Microchiroptera) that inhabit the cave, littered the surface of the 2004 tsunami deposit. The basal contact of the 2004 deposit is an erosional unconformity.
Prehistoric tsunami deposits. Beneath the 2004 tsunami deposit, we found an additional 11 sand beds (A-K) that we interpret as tsunami deposits (Fig. 4). The 11 sand beds consist of well-sorted, normally graded, very fine sand to silt with a sharp basal contact. There is no evidence of unconformities in the stratigraphic sequences from the trench-wall exposures ( Supplementary Figs 1 and 2). Sand beds G-J have thin deposits (2-7 cm), whereas sand bed F has the thickest deposit (23 cm). Some of the sand beds have a rip-up clast-rich lower portion and a lenticular-laminated upper portion. The rip-up clasts are very similar to the deposits that underlie them. Large detrital weathered fragments of cave chalk are preserved in the sand beds. Foraminifera are abundant, in particular in sand beds I-K (Fig. 5;  Supplementary Tables 3 and 4). The provenance of the foraminifera ranges from intertidal to subtidal to offshore 23 . A large percentage of the foraminiferal assemblage in each sand bed is pristine (Fig. 5).   Intercalated beds. Organic-rich beds are found between the 11 sand beds (A-K), reflecting slow sediment accumulation during intervals between tsunamis. The organic beds are commonly finely laminated and range in thickness from o1 mm to 9 cm ( Fig. 4; Supplementary Figs 1 and 2). The intercalated beds consist of sands, reworked by periodic drips of water through the cave ceiling during periods of high precipitation and insect burrowing. The organics were likely produced by the same processes that produce organic debris on the surface of the 2004 tsunami and have been broken down by post-depositional processes. Foraminifera are absent or in low abundances with a fragmented and abraded assemblage of intertidal to subtidal to offshore species (Fig. 5), further suggesting the intercalated beds are reworked from the tsunami sand beds A-K. In many intercalated beds, we found small, pristine and fragile chalk florets.
Four mud beds appear between sand beds B-H with thicknesses up to 25 cm (Figs 4 and 5; Supplementary Figs 1 Tables 3 and 4). However, foraminifera are absent in the other two mud beds, suggesting deposition by freshwater ponding in the cave. Pristine cave chalk is found in some mud beds, further supporting deposition in low-energy conditions.
Chronological constraints. We obtained accelerator mass spectrometry (AMS) radiocarbon ages on pieces of detrital charcoal and whole molluscs from within, below and above the sand beds ( Fig. 6; Table 1). We interpret the two bracketing dates as maximum and minimum ages for the timing of sand bed deposition. Fragments of charcoal from an organic-rich bed at the base of the sedimentary sequence yield a maximum age of 7,672-7,588 years BP for sand bed A. Charcoal from a mud bed (5,583-5,

Discussion
Coastal caves have not previously yielded prehistoric records of tsunamis. Indeed, the cave's sheltered location and absence of human activity suggest that these sand beds represent the best-preserved and most complete tsunami history for the Indian Ocean between 7,400 and 2,900 years BP. The cave's interior protects the tsunami deposits from erosion. The rock sill near the cave entrance (Fig. 2) mitigated the erosional impact of tsunamis that are found at elevations beneath the sill. However, deposits above the rock sill are vulnerable to scouring from subsequent events. The cave's location also disfavours sand bed deposition or re-working by intense storms 12,14 . Exposure to tropical cyclones is limited due to the lack of Coriolis force near the equator 24,25 .
In addition, the track of any tropical cyclones that originate in Indian Ocean will move towards India, Bangladesh or Myanmar without producing a storm surge in Sumatra 14,26 . Although tropical cyclones do strike eastern Thailand, they dissipate after crossing the Malay Peninsula and Sumatra before moving offshore along Sumatra's west coast (for example, tropical storm Vamei in 2001 (refs 27,28)). The stratigraphic and microfossil data of the 11 prehistoric sand beds (A-K) resemble the 2004 tsunami as well as tsunami deposits described elsewhere. Rip-up clasts at the base of sand beds and sharp basal contacts suggest erosion occurred at the beginning of the tsunami inundation as the surge entered the cave. Normally graded sand beds indicate settling from suspension following tsunami inundation in the cave 9,29 . The normal grading suggests that each bed resulted from a single (rather than multiple) instance of the cave filling with water and draining. The foraminifera assemblage of the sand beds were dominated by intertidal to subtidal to offshore species. Marine foraminifera often dominate tsunami deposits because of the landward transport and deposition of scoured marine sediment 23,30 . The taphonomic (or surface) condition of individual foraminifera distinguishes the tsunami sand beds and the intercalated beds (Fig. 5). The foraminifera of the tsunami sand beds is predominantly pristine suggesting the foraminifera were entrained from a protected subtidal substrate 31,32 .
We have also identified cave chalk weathering as a new indicator of tsunami inundation. Large fragments of weathered cave chalk are preserved in the sand beds. These fragments most likely fell from the cave ceiling and were weathered due to abrasion by tsunami transport. In contrast, we found pristine and fragile chalk florets in the organic and mud beds.
Holocene relative sea-level reconstructions from the Indo-Pacific region are characterized by a mid-Holocene sea-level high stand of  a few decimetres to several metres 33,34 , but the presence or absence of such a highstand may be controlled by local tectonic processes 35 . A record of buried soils from the northwestern coast of Aceh Province suggest that relative sea-level rose during the early and mid-Holocene from À 5 m at B7,900 years BP to À 1.6 m at B5,700 years BP 12 . Relative sea-level was below present until at least 3,800 years BP. In the late Holocene, relative sea-level stabilized within 0.4 m of modern sea-level 12,22,33 . This gradual long-term relative sea-level rise without a mid-Holocene highstand created a time-window for tsunami deposits and intercalated beds to aggrade without a significant interruption in sedimentation 36 (Fig. 7). The cave probably contained stratigraphic evidence of recent historic tsunamis from 2,900 years BP to the 2004 Indian Ocean tsunami that have been identified elsewhere in the region 10,14 , but these were most likely removed by subsequent tsunamis inundating the cave as indicated by the erosional unconformity beneath the 2004 deposit (Figs 4 and 7; Supplementary Fig. 1).  The missing stratigraphic record coincides with the continued aggradation of the nearby coastal plain 12 . Stratigraphical records of late Holocene tsunamis are generally restricted to environments with sufficient accommodation space, such as intervening coastal swales between ridges 10,14 , estuaries and ponds where overwash deposits are protected from erosion by rapid growth of vegetation or deposition of sediment 36 .
Independent evidence for tsunami inundation in the cave comes from stratigraphy of nearby coastal lowlands 12,37 . Three coseismic subsidence events and seven tsunamis between B7,500 and 3,800 years BP are documented in the stratigraphy of the west coast of northern Sumatra 12 . During this time interval, the cave sequence preserves an identical number of tsunamis (that is, sand beds A-G; Fig. 6 and Table 1). Offshore of Sumatra, Patton et al. 38 identified 11 deep-sea turbidites along the Andaman-Aceh slip patch between B6,500 and 2,700 years BP. Although the number of events is the same, the timing of the events is different. The deep-sea turbidite record does not capture the tightly clustered tsunamis (sand beds G-J) and the large gap in time between tsunamis F and G. The discrepancies suggest that ruptures along the Sunda megathrust do not always trigger both tsunami deposits and turbidites.
Since the 2004 tsunami, considerable evidence for prehistoric tsunamis has been obtained from sites around the Indian Ocean [13][14][15][16]18,20,21,[39][40][41][42][43] . However, studies with time spans comparable to the cave are restricted to Sri Lanka 19 and the Maldives 21 . In southern Sri Lanka, Jackson et al. 19 , identified seven tsunami sand beds between B6,700 and 2,400 year BP. Klostermann et al. 21 identified three tsunami sand beds between B5,600 and B2,900 years BP in the Maldives. However, these far-field records do not capture the tightly clustered tsunamis (sand beds G-J) and have events that span the large gap in time between tsunami sand beds F and G. The far-field seismic sources for these tsunamis are uncertain. For example, slip along the megathrust near the Andaman-Nicobar Islands are potential seismic sources for tsunamis in Sri Lanka and the Maldives 44 . In addition, the faults along the southern coast of Pakistan are potential seismic sources for tsunamis in the Maldives 44 . The immediate proximity of the cave to the Sunda megathrust provides a more reliable indicator of tsunamis generated by ruptures of the megathrust than far-field records.
The chronology from accelerator mass spectrometry (AMS) radiocarbon ages from the sand beds of the coastal cave and the stratigraphy of the nearby coastal lowlands 12,37 , combined within a  Table 5).
Bayesian framework ( Fig. 6; Methods section), provides the chronology of tsunamis between 7,400 and 2,900 years BP ( Table 1). The chronology suggests an average recurrence interval of 456 years between 7,400 and 2,900 years BP with a large uncertainty (95% C.I. 1-2,210) ( Fig. 6; Supplementary Table 7). A similar average recurrence interval (600-900 years) was estimated from the nearby coastal lowlands of northwestern Aceh Province between 7,400 and 3,800 years BP 12 . The tsunami record from Sri Lanka 19 suggested an average recurrence interval of B360 years between 6,600 and 4,200 year BP. The tsunami record from the cave, however, indicates a dramatic variation in recurrence interval. Between 7,400 and 5,500 years BP, the recurrence interval for tsunamis A to F was 681 years (95% C.I. 11-2,222) (Fig. 6). But after 5,500 years BP, the coastal cave has an age gap of 2,164 years (95% C.I. 1,997-2,247) between tsunamis F and G. Four tightly clustered tsunamis (G-J) occurred between 3,400 and 3,300 years BP with an average recurrence interval of 16 years (95% C.I. 0-55). The most recent tsunami (K) recorded in the coastal cave occurred at 2,900 years BP, with a recurrence interval of 426 years (95% C.I. 357-505). Although the time span of the northern Simeulue coral microatoll record is limited to the last millennium, it shows a similar large variation of recurrence intervals from 56 years to B550 years 39 .
There is a correlation between the thickness of tsunami sand beds and recurrence intervals in the cave (Fig. 6). The thinner sand beds (G-J) have the smallest recurrence intervals and were preceded by the largest age gap between tsunamis (Supplementary Fig. 4). The thickest sand bed (F) preceded the large age gap, with a thickness similar to the 2004 tsunami sand bed. Although variations of offshore sediment availability or lateral shoreline changes might play a role in the thickness of tsunami beds 45 , we suggest that the thickness of the sand beds may reflect the size of slip along the megathrust. It is possible that sand bed F was deposited by a giant tsunami produced by a large slip that was followed by a very long dormant, interseimic period with substantial strain accumulation. Subsequently, partial, smaller slip failures occurred in rapid succession between 3,400 and 3,300 years BP, producing sand beds G-J. The very long dormant period suggests that the Sunda megathrust is capable of accumulating large slip deficits between earthquakes. Such a high slip rupture would produce a substantially larger earthquake than the 2004 event.
The dramatic variation in tsunami recurrence intervals suggests a continuum of recurrence behaviour from large slip ruptures (2004 tsunami and sand bed F), earthquake super cycles or doublet earthquake 2,4,11,39 to smaller slip failures (for example, sand beds G-J) similar to the October 2010 Mentawai tsunamigenic earthquake (Fig. 1). Variations in recurrence may result from temporal changes in coupling or locking depth, or very long-term slow non-tsunamigenic slip events 46,47 . If thickness of the tsunami deposit ( Fig. 6) reflects slip and the size of the slip patch, the thickness of the 2004 tsunami deposit implies a long dormant period until the next large slip event. Slow non-tsunamigenic slip events might predominate during such long periods of quiescence and precede clustering smaller slip failures along the megathrust. The remarkable variability of recurrence suggests that regional hazard mitigation plans should be based upon the high likelihood of future destructive tsunami demonstrated by the cave record and other paleotsunami sites, rather than estimates of recurrence intervals.

Methods
Field methods and data collection. The evidence for prehistoric tsunamis is derived principally from stratigraphic relations found in six trench-wall exposures and other smaller pits in the cave. The stratigraphy is visually striking, because of a strong contrast in colour between alternating beds of sand, mud and laminated organic sand ( Supplementary Figs 1 and 2).
We identified and described stratigraphic units in Trench 1 and the trench walls were mapped from high-resolution photos of the vertical surfaces. The alternating beds of sand, clay and laminated organic sand allow stratigraphic correlation of individual beds within the trench-wall exposure ( Supplementary Figs 1 and 2). We divided the stratigraphy into 12 units, from Unit 1, which is the oldest and deepest,  ARTICLE resting on the limestone bedrock of the cave, through Unit 12, which is the youngest unit that underlies the 2004 Indian Ocean tsunami sand bed (Fig. 5).
Column samples were taken across each stratigraphic interface, providing samples for grain size, microfossil and chronology. Topographic data were based on a total station survey in 2011 and 2012 (Fig. 2). Excavation of six trenches within the cave interior was conducted in 2011 and 2012. We made measurements of sedimentology and stratigraphy in the field.
Tsunami sand beds. We use several lines of stratigraphic evidence that point to a tsunami origin for sand beds A-K that includes stratigraphic relations, lithology, degree of sorting, internal structures, such as fining upwards beds, rip-up clasts of mud, sharpness of stratigraphic contacts, uniformity in bed thickness and lateral continuity of beds 8,[48][49][50] . In addition, we use secondary features such as liquefaction, bedding-plane faults (décollements) and normal faults that are perhaps triggered by seismic shaking 8 .
We analysed foraminiferal assemblages ( Fig. 5; Supplementary Tables 3 and 4) from most units to confirm marine inundation and indicate the provenance of the overwash sands 32 . In addition, we used foraminiferal test taphonomy 51 to reveal depositional and post-depositional environmental conditions (Supplementary Fig. 3; Supplementary Tables 3 and 4). The foraminiferal analyses help constrain sediment provenance and help identify tsunami sands in a variety of environmental settings.
We used the taphonomy of cave chalk as a new indicator of the environment of deposition. The cave chalk found in the sedimentary units most likely fell from the cave ceiling due to weathering. The chalk occurs in two forms: (1) small, fragile, pristine florets; and (2) large detrital rounded fragments that do not have a delicate floret structure (Supplementary Fig. 3). Since the florets do not show evidence of water transport (for example, rounding), we suggest that the delicate florets are limited to sand beds not deposited by tsunamis. In contrast, the rounded chalk fragments are only found in sand beds that were produced by a tsunami. We believe the rounding of the fragments is due to abrasion by the tsunami waves.
Grain size and thickness data allows comparison of the 2004 tsunami sediments with the prehistoric sand beds (Figs 5 and 6; Supplementary Table 1). Further, the general fining upwards trends in many tsunami beds suggests deposition by tsunami waves (Fig. 5; Supplementary Table 1). Following an initial high water flow, a decreased flow velocity often causes sand beds to deposit in graded, fining upwards sequences 9,52 . For example, along the Aceh coastline, the 2004 tsunami sand bed records distinct fining upwards sequences of sand deposition 52 , similar to the 2004 tsunami sand bed preserved in the coastal cave. Samples for foraminifera, cave chalk and grain size analyses are from Trenches 1 and 4 ( Fig. 5; Supplementary Tables 1, 2, 3 and 4).
Radiocarbon dating. We collected detrital charcoal fragments from seven stratigraphic units (Units 1, 5, 8, 9, 10 and 11), and two intact gastropods from Unit 6. The radiocarbon ages constrain ages of tsunamis; calibrated radiocarbon ages helped constrain the timing of tsunamis. Detrital charcoal and gastropods were collected from units for radiocarbon dating and were analysed by GNS, Rafter Radiocarbon Laboratory, New Zealand (Table 1).
We calibrated radiocarbon ages with Calib rev. 6.0.0 (ref. 53). The calibrated age ranges appear with 95.4% HDR (B2 standard deviations), where years 'before present' (BP) is years before A.D. 1950 (Table 1). Also, we corrected the ages of the gastropod shells for the marine reservoir effect using a DR value of 15 ± 119 (ref. 54) to account for the fact that Indian Ocean waters show substantial 14 C depletion due to upwelling.
To further constrain chronology, we analysed the organic guano beds derived from insect-feeding bats that occupy the coastal cave. We sampled six organic-rich beds along thin stratigraphic horizons (3-6 mm). Although the radiocarbon analyses indicate the dark sand beds are broadly mid-late Holocene in age, three radiocarbon dates are not consistent with stratigraphic position. We suggest that the discrepancies of the bulk guano dates are due to: (1) bulk samples containing an unknown mixture of organic material of variable age, and representing an average age of the sample; and (2) groundwater percolating along cracks in the limestone cave, and introducing exogenous, old carbon into the organic beds.
Bayesian age-depth model. We apply a Bayesian age-depth model to 19 radiocarbon dates from the coastal cave and the nearby coastal lowlands 12 . The radiocarbon dates either directly date a tsunami or provide maximum or minimum age limits for a tsunamis ( Fig. 6; Table 1) (Supplementary Tables 6  and 7). Our Bayesian modelling approach provides control over the model fitting process and flexibility in the modelling assumptions. The code is available at https://github.com/andrewcparnell/tsunamis. We use the following notation to build our model. y i is the calendar age of tsunami i, where i runs from 1 to 11. These are the parameters we are most interested in estimating. Together, we write these values as y. x ij is the direct radiocarbon date j of tsunami i, where j ¼ 1,y,n i with n i the number of direct dates for tsunami i. These values have associated fixed 1-sigma errors s ij . Note, that for some tsunamis there are no direct dates, in which case n i ¼ 0. Thus, while we have 13 direct dates in total, five tsunamis are without direct dates. Together, we write these values as x. y Ã i are the calendar ages of limiting dates lying between tsunamis i and i þ 1. They are nuisance parameters. Together, we write these values as y Ã . x Ã ij is a limiting radiocarbon date j lying between tsunamis i and i þ 1. These values have associated fixed 1-sigma errors s Ã ij . Here j ¼ 1; . . . ; n Ã i , where n Ã i represents the number of limiting dates for tsunami i. As above, some of these are 0. Together, we write these values as x Ã . g i are the calendar age shifts of the limiting dates which provide a maximum age of tsunami i. They are nuisance parameters which we write together as g. x ÃÃ ij is a limiting radiocarbon date j for tsunami i, providing evidence of a maximum age. These values have associated fixed 1-sigma errors s ÃÃ ij . Here j ¼ 1; . . . ; n ÃÃ i ; where n ÃÃ i represents the number of limiting dates for tsunami i. Again, some of the n ÃÃ i values are 0, as we have only three such dates. Together we write these values as x ÃÃ . Note that the calibrated value of a radiocarbon age x ÃÃ ij is y i þ g i , that is, the calendar age of the tsunami plus a shift indicating how much older the radiocarbon date is beyond that of the tsunami itself. r(y) is the IntCal13 calibration curve which has the probability distribution r(y)BN(m(y),t 2 (y)).We assume that both m() and t 2 () are known functions.
Our overall goal is to find the posterior distribution: where are the likelihood terms. The prior distribution p(y) is set to enforce the ordering of the dates: where I is an indicator function. The prior distribution on the excesses g i is the only informative prior in the model. We use: which corresponds to a gamma distribution with a mean of 200 and a s.d. of 200. While this distribution is diffuse, it is informative for the maximum ages. The distribution corresponds to a 95% probability that the maximum ages are no more than 600 years older than the tsunami they are aiming to represent. Our model is fitted with Markov chain Monte Carlo using Metropolis-Hastings steps for all parameters since ordering the tsunami dates was complicated by their prior distribution 55 . The model is sensitive to starting values of the parameters due to the ordering constraint, so we simulate suitable values by calibrating dates individually, and sampling from these distributions with an extra restriction on the ordering. We run the model with multiple different starting values, and check convergence using trace plots and the Geweke convergence diagnostic 55 . The final model run created 1 million iterations, removing 100,000 for a burn-in period and then keeping only every 450th iteration.
2004 tsunami sand bed. The 2004 tsunami sand bed ranges from about 20 to 43 cm in thickness in trench-wall exposures ( Fig. 3; Supplementary Figs 1 and 2). It is a light grey, normally graded fine to very fine sand (mean ¼ 2.7F; % sand ¼ 91.9%), with abundant laminations, some of which can be traced more than a metre. In Trench 6, the 2004 tsunami sand bed records three distinct beds (Fig. 3), delineated by three pulses of coarse material followed by subsequent fining upwards sequences (Fig. 3; Supplementary Table 1). We interpret each combined coarse pulse and fining upwards sequence as individual tsunami waves. The first coarse pulse at the base of the 2004 sand (41-42 cm) is marked by an influx of fine to medium sand (mean ¼ 2.1F; % sand ¼ 94.1%), which fines upwards until 34 cm (mean ¼ 2.7F; % sand ¼ 90.6%). The second coarse pulse occurs from 31-33 cm (mean ¼ 2.3F; % sand ¼ 95.6%) and fines slightly up to 24 cm (mean ¼ 2.7F; % sand ¼ 92.5%). The final coarse pulse is located between 22 and 24 cm (mean ¼ 2.5F; % sand ¼ 96.3%), and fines upwards to a very fine sand (mean ¼ 3.7F; % sand ¼ 74.0%) at the top of the sequence (0-1 cm).
Rip-up clasts, consisting of organic-rich granules, wood and shells are common, especially in the lower part of the 2004 deposit. Abundant (2,500-3,246 individuals per 1 cm 3 ) foraminifera that are predominantly pristine (41-52%) and sourced from subtidal (52-58%), intertidal (33-38%) and offshore (for example, planktic) (8-14%) environments are present, as are weathered fragments of cave chalk (Fig. 5, Supplementary Fig. 3; Supplementary Table 4). The 2004 sand has the highest diversity foraminiferal assemblage, with Pararotalia sp., Amphistegina sp. and Calcarina sp. dominating. A sharp and erosional contact marks the boundary between the 2004 tsunami sand and the underlying Unit 11. The erosional removal of pre-2004 sediment is variably preserved in trenches and on the cave walls as alternating remnants of sand and organic-rich sand.
Unit 1 sand bed A. The lowest stratigraphic unit (Unit 1) above the limestone cave floor is an irregular, laminated, dark, organic layer. Overlying the Unit 1 organic layer is fine sand (Unit 1.1; mean ¼ 3.0F; % sand ¼ 77.4%) that does not contain foraminifera or cave chalk. Unit 1.2 is discontinuous marine-influenced clay (grain size data not available) that pinches out toward the corner of Trench 1. The clay bed contains relatively low numbers of foraminifera (20 individuals per 1 cm 3 ). The species assemblage is dominated by intertidal (71%) and subtidal (29%) species, with a paucity of planktic and deeper-dwelling benthic foraminifera. The taphonomic assemblage is dominated by abraded (52%) and fragmented (48%) individuals, with no pristine individuals present ( Fig. 5; Supplementary Table 1). Unit 1.2 is devoid of cave chalk. Unit 1.2 is overlain by a 5.4 cm thick fine sand (Unit 1.3; mean ¼ 2.9F; % sand ¼ 74.9%) with sparse laminations and abundant rip-up clasts derived from the underlying clay. The lower stratigraphic contact between Unit 1.3 and the underlying clay is sharp (B2 mm), along an erosional and irregular surface. On the basis of the abundance of rip-up clasts, we interpret Unit 1.3 as the oldest tsunami bed, labelled Sand Bed A. Unit 1.3 sand contains a low number of foraminifera (21-62 individuals per 1 cm 3 ) consisting predominantly of subtidal (38-58%) and intertidal (41-62%) species. The assemblage is dominated by the subtidal species Epinoides sp., Cibicides lobatulus and Pararotalia spp., with a near absence of planktics. Individual foraminifera were both abraded (68-77%) and fragmented (14-30%), with only 2-9% in pristine condition ( Fig. 5; Supplementary Table 3  Overlying Unit 2.1 is a B6 cm thick fining upwards, fine to very fine sand (Unit 2.2; mean ¼ 3.5F; % sand ¼ 75.2%; Supplementary Fig. 1; Supplementary  Table 1). Along the base of Unit 2.2, pebble-sized clay rip-up clasts are abundant. Locally, Unit 2.2 grades into a thin, discontinuous dark sand. Unit 2.2 contains weathered cave chalk and moderate abundances of foraminifera (76-89 individuals per 1 cm 3 ) from subtidal (41-47%) and intertidal (51-59%) environments ( Fig. 5; Supplementary Tables 3 and 4). Unlike the underlying Unit 2.1 foraminiferal-bearing sediments, species within Unit 2.2 (Fig. 5) include the deeper-dwelling Lagena sp., and Operculina ammonoides. The rip-up clasts, normal grading and the presence of subtidal foraminifera and weathered cave chalk in Unit 2.2 imply deposition by a tsunami (sand Bed B). However, the foraminiferal assemblage of sand bed B, unlike most other candidate tsunami sand beds in the cave sequence, is mostly fragmented (47-55%), and contains only minor abundances of pristine individuals (9-12%). The contact between Unit 2.2 and the overlying Unit 3 clay is sharp (B2 mm), but irregular.
Unit 4 sand bed D. The base of Unit 4 consists of a 4-cm thick, laminated, organic-rich sand. The laminations are clear but discontinuous, and vary from black to dark red-brown. Thickness also varies markedly, from o1-12 cm. The sand lacks foraminifera and cave chalk. Overlying Unit 4 is a 3-cm thick, fine sand (Unit 4.1; mean ¼ 2.5F; % sand ¼ 86.8%; Fig. 5, Supplementary Table 1). The stratigraphic contact between Unit 4.1 fine sand and the underlying Unit 4 is sharp (a few mm), along an irregular surface. Unit 4.1 contains a low number of foraminifera (19 individuals per 1 cm 3 ). Many of the tests are abraded (78%) and fragmented (22%), and no pristine individuals were found. Where present, foraminifera were of intertidal (83%; Ammonia convexa, Ammonia parkinsoniana) or subtidal (17%; Pararotalia sp.) origin ( Fig. 5; Supplementary Table 4). No cave chalk was observed in Unit 4.1. On the basis of the presence of fine sand, with a basal erosional contact and subtidal foraminifera in Unit 4.1, we infer a speculative tsunami, labelled Sand Bed D (Supplementary Fig. 1). Overlying Unit 4.1 is a clayey mud (Unit 5) with a sharp (B2 mm), but irregular contact.
Unit 5 sand bed E. Unit 5 consists of an extensive non-marine clayey mud. The clayey mud is 11 cm thick, but varies in thickness from 4 to 15 cm. The massive clay contains pristine cave chalk, but no foraminifera and was likely deposited by ponding water in the cave resulting from periods of increased precipitation. The upper contact of the clay is smooth and displays only minor, local evidence of erosion. Overlying the clayey mud is a massive 3-cm thick, very fine sand (Unit 5.1; mean ¼ 3.4F; % sand ¼ 72.9%), that has a low foraminiferal assemblage (35-42 individuals per 1 cm 3 ; Fig. 5, Supplementary Tables 1 and 3). The foraminifera are highly abraded (63-66%), and predominantly from an intertidal (Ammonia parkinsoniana and Elphidium advenum) environment (56-62%). Unit 5.1 is devoid of cave chalk. On the basis of the presence of massive sand and marine foraminifera, we speculate that Unit 5.1 (labelled Sand Bed E) represents deposition as a result of a tsunami. The absence of rip-up clasts from the underlying clay suggests that the flow of the water in the cave was too weak to erode the underlying clay. In addition, it is hard to explain how the massive sand could have been deposited from other processes so far back into the cave. Overlying Unit 5.1 is a 3-cm thick, laminated, very fine sand that varies from black, guano-rich lamina to reddish brown, inorganic lamina (Unit 5.2: mean ¼ 4.0F; % sand ¼ 58.3%). Along other trench wall exposures, a thin, massive, fine sand (Unit 5.3: no grain size data available) overlies Unit 5.1. The fine sand lacks rip-up clasts, and contains fragments of pristine cave chalk and a very low number (2 individuals per 1 cm 3 ) of abraded (100%) foraminifera, representing reworking of older units. The contact between Unit 5.1 (sand bed E) and the overlying units is irregular. Overlying the thin organic-rich fine sand (Unit 5.2) are the thickest non-marine clayey mud beds exposed in the cave sedimentary sequence (Units 5.4 and 5.5; no grain size data available) ( Supplementary Figs 1 and 2). Their combined thickness ranges up to about 25 cm. The clayey mud is differentiated into two separate units based upon their colour: the upper part of the clayey mud (that is, Unit 5.3) is darker (more organic) than the lower part of the clay (Unit 5.4). Both units are devoid of foraminifera and were likely deposited by the same processes that formed Unit 5. The units lack cave chalk.
Fragments of detrital charcoal from the overlying Unit 5.4 clayey mud yielded a calibrated age range of 5,583-5,331 years BP, which we interpret as the minimum age of sand bed E.
Further evidence from secondary folding of Units 5.1 and 5.2 and the truncation of these units along the top of Unit 5.3 clayey mud suggests local ground shaking ( Supplementary Fig. 2c). On the basis of the intact stratigraphic continuity of Units 5.3 and 5.4, ground shaking must have occurred before the deposition of Unit 5.3.  Tables 3 and 4), with Elphidium advenum, Calcarina spp. and Operculina ammonoides dominating the assemblage. While the Unit 6 lower sand contains weathered fragments of cave chalk, the upper sand is devoid of them. The abundance of rip-up clasts, graded bedding, subtidal foraminifera and weathered cave chalk suggests a tsunami origin for Unit 6, labelled sand bed F. Similarly to other sand units that represent tsunami deposits, including the 2004 tsunami deposit, the majority of foraminifera are generally pristine (41-69%). Overlying Unit 6 ( Supplementary Figs 1 and 2; sand bed F) is a o2.0-cm thick, black, organic sand that grades into a heterogeneous, laminated, grey sand (Unit 7). The contact between Unit 6 and the overlying Unit 7 is sharp (a few millimetres) and irregular.
Two pristine gastropod shells from Unit 6 yielded a calibrated ages range of 5,231-4,515 years BP and 5,258-4,552 years BP, which we interpret as the age of Sand Bed F. Unit 7 sand bed G. Underlying Unit 7.1 or Sand Bed G is a 1.5-cm thick, black, organic, fine sand (Unit 7: mean ¼ 2.6F; % sand ¼ 89.4%) with discontinuous laminations. Its thickness varies from a few millimetres to about 5 cm. The contact between Unit 7 and the overlying Unit 7.1 is sharp (B2 mm) along an erosional surface. Unit 7.1 consists of irregularly laminated fine sand (mean ¼ 3.0F; % sand ¼ 79.8%) that ranges from 1 to 4 cm in thickness. Unit 7.1 fines upwards from fine sand at the base (mean ¼ 2.6F; % sand ¼ 85.1%) to very fine sand at the top (mean ¼ 3.6F; % sand ¼ 72.8%; Supplementary Fig. 2; Supplementary  Tables 1 and 2). Unit 7.1 contains weathered cave chalk and foraminifera (93-105 individuals per 1 cm 3 ) that are dominantly pristine (65-82%) and subtidal (71-77%; Supplementary Fig. 2; Supplementary Table 4). Pararotalia sp. and Cibicides spp. dominate the foraminiferal assemblage. In places, the fine sand laminations are defined by heavy minerals. The normal grading, the presence of subtidal foraminifera and weathered cave chalk imply deposition by a tsunami (Fig. 5 sand bed G). Sand bed G (Unit 7.1) is overlain by a thin clayey mud (Unit 7.2: mean ¼ 4.5F; % clay ¼ 14.7; % silt ¼ 34.6) that is laterally discontinuous. Unit 7.2 is devoid of both foraminifera and cave chalk fragments. The contact is gradational over a few millimetres. The abundance of graded bedding, subtidal foraminifera and weathered cave chalk suggests a tsunami origin for Unit 8, labelled sand bed H ( Supplementary  Figs 1 and 2). Overlying Unit 8 or Sand bed H is a thin, clayey mud bed (Unit 8.1: mean ¼ 4.6F; % sand ¼ 53.5%) that is only exposed in Trench 4. The contact between Unit 8.1 and the underlying sand bed H (Unit 8.1) is irregular, but sharp over a few millimetres. Unit 8.2 consists of massive to laminated, very fine sand (mean ¼ 3.2F; % sand ¼ 82.9%). Unit 8.2 is devoid of foraminifera and cave chalk. Unit 8.3 consists of thin, black and grey to dark brown, organic, sandy-mud (mean ¼ 5.7F; % sand ¼ 19.8%), about 2 cm thick. The dark organic sandy-mud contains pristine cave chalk fragments and a low number (20 individuals per 1 cm 3 ) of foraminifera that are generally abraded (96%), and sourced from intertidal (55%) and subtidal (45%) environments. The contact between Unit 8.2 and the overlying Unit 9 sand is sharp and irregular.
Further evidence from a secondary high-angle, normal fault that disrupts Units 1 through 7 ( Supplementary Fig. 1e), suggest that faulting occurred before deposition of Units 8 or sand bed H.
Fragments of detrital charcoal from Unit 11 yielded a calibrated age range of 3,362-3,246 years BP, which we interpret as the age of sand bed H. Unit 9 sand bed I. Unit 9, consists of a light tan, massive, fine sand (mean ¼ 2.3F; % sand ¼ 86.3%), up to 3 cm thick, with weathered detrital fragments of cave chalk, particularly in the upper half of the unit (Supplementary Fig. 2; Supplementary Table 4). Angular fragments of detrital charcoal are also common. Faint laminations are present, in particular at the base of the unit. Unit 9 does, however, locally fill small scours cut into the underlying Unit 8.3 sandy-mud. This indicates minor erosion prior to deposition. Foraminifera are abundant (256 individuals per 1 cm 3 ), and are predominantly subtidal (58%) species (Pararotalia stellata, Epinoides repandus and Cibicides refulgens); B44% of the foraminiferal assemblage is pristine ( Fig. 5; Supplementary Tables 3 and 4). The presence of fine sand, subtidal foraminifera and weathered cave chalk imply deposition by a tsunami, labelled sand bed I (Fig. 5 band bed I). Unit 9 sand (band bed I) is overlain by a black, organic bed only a few millimetre thick (Unit 9.1: no grain size data), with a sharp contact. A low-abundance assemblage (4 individuals per 1 cm 3 ) consisting of abraded (82%), intertidal (54%), and subtidal (46%) foraminifera and fragments of pristine cave chalk is present in this unit, and were derived from underlying band bed I (Unit 9). Fragments of detrital charcoal from Unit 9 yielded two calibrated age ranges of 3,358-3,221 years BP and 3,366-3,252 years BP, which we interpret as the age of sand Bed I. In addition, fragments of detrital charcoal from the top of Unit 9 sand yielded a calibrated age range of 3,363-3,245 years BP, which we interpret as the maximum age of Sand Bed I.
Further evidence from truncation of Units 8.1, 8.2 and 8.3 suggests a bedding-parallel fault along the top of Unit 8 ( Supplementary Fig. 1c). Secondary folding of these units above this bedding-parallel fault implies that displacement occurred after deposition of Unit 8.3 but before deposition of Unit 9 (sand bed I).
Unit 10 sand bed J. Unit 10 is a thin, very fine sand with red-brown laminations (mean ¼ 3.2F; % sand ¼ 75.6%). Unit 10 contains a low abundance (nine individuals per 1 cm 3 ) and weathered (76% fragmented and 24% abraded) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms16019 foraminiferal assemblage (86% intertidal species), with pristine cave chalk fragments. This unit represents accumulation of bat guano and remobilized sand from the underlying Unit 9 sand. Overlying the red, very fine sand of Unit 10, along a sharp contact, is a thin (2 cm thick), massive fine sand (Unit 10.1; mean ¼ 2.8F; % sand ¼ 85.9%; Supplementary Figs 1 Tables 1 and 2), with abundant pebble-sized fragments of weathered cave chalk. Unit 10.1 sand contains abundant (190 individuals per 1 cm 3 ) foraminifera that are predominantly pristine (53%) and from subtidal (51%) and intertidal (45%) environments, with only a small percentage of planktic species present (4%; Supplementary Fig. 2; Supplementary Tables 3 and 4). Dominant species include Ammonia parkinsoniana, Elphidium craticulatum and Heterolepa sp. The contact with the underlying Unit 10 is sharp. On the basis of the presence of very fine well-sorted sand and subtidal foraminifera, we infer that Unit 10.1 ( Fig. 5; labelled sand bed J) represents deposition as a result of a tsunami. Overlying Unit 10.1 (sand bed J) is a black to dark brown, massive, organic, fine sand (Unit 10.2; mean ¼ 2.5F; % sand ¼ 87.8%). Unit 11 has been eroded away from most of the other trench exposures. Unit 10.2 is thin, black to dark brown, massive, organic, fine sand, about 1.5 cm thick, that represents bat guano accumulation and sand remobilization derived from Unit 10.1.

and 2; Supplementary
Fragments of detrital charcoal from Unit 10.1 yielded four calibrated ages: 3,464-3,383 years BP, 3,370-3,260 years BP, 3,356-3,219 years BP, 3,236-3,068 years BP. These calibrated ages represent the age of sand bed J. In addition, fragments of detrital charcoal from Unit 10 yielded a calibrated age range of 3,396-3,269 years BP, which we interpret as the maximum age of sand bed J.
Unit 11 sand bed K. Unit 11 is an 8.4-cm thick sequence that fines upwards from a medium sand (Unit 11; mean ¼ 1.9F; % sand ¼ 95.5%) at the base to a fine sand (mean ¼ 2.4F; % sand ¼ 94.3%) with numerous laminations at the top (Supplementary Figs 1 and 2; Supplementary Tables 1 and 2). The unit contains abundant rip-up clasts that are derived from the underlying Unit 10.2. The uppermost laminated, very fine sand is nearly absent in other excavations due to erosion, but it is preserved in most of the Trench 1 exposures ( Supplementary  Fig. 1). Weathered fragments of cave chalk are abundant and range in size up to 3 cm. Unit 11 also abounds in foraminifera (1,976-2,485 individuals per cm 3 ), and 49-55% of the foraminifera are pristine. Subtidal (48-61%) species are dominant (Pararotalia sp., Asterorotalia sp. and Epinoides repandus) with planktic foraminifera present, but in smaller amounts (9-14%; Fig. 5, Supplementary  Tables 3 and 4). On the basis of the normal grading, abundant rip-up clasts and subtidal foraminifera, we suggest that this sand layer represents a tsunami deposit ( Fig. 4; labelled sand bed K). Sand bed K (Unit 11) is separated from the overlying black to grey, thin, laminated, fine sand (Unit 12: no grain size data available) by an erosional unconformity. Unit 12 is a massive, fine sandy pebble breccia with an irregular black, organic sand cap. Clasts range up to about 10 cm in length and are commonly angular.
Fragments of detrital charcoal from Unit 11 yielded two calibrated age ranges of 2,859-2,772 years BP and 2,975-2,862 years BP. We interpret these as the age of sand bed K.
Grain size analysis. We analysed grain size at 1 cm resolution with a Malvern MS 3000 laser particle size analyzer (measuring grain sizes up to 1,800 mm). Before analysis, we removed organics and carbonate with 30% hydrogen peroxide and 10% hydrochloric acid, respectively. We subsequently let the samples disaggregate in a sodium hexametaphosphate solution for 24 h.
We calculated grain size values with the Wentworth-Phi scale 56 , using the average of three runs. Grain size descriptions for each sampled interval follow those defined by Blott and Pye 57 , and include a mean (average grain size), mode (dominant grain size), s.d. (degree of sorting) and percentage of clay, silt and sand. We show the depths of coarse pulses and fining upwards sequences in Fig. 5.
Foraminiferal and cave chalk analysis. We examined 42 samples from Units 1-12 for foraminiferal taxa (Supplementary Table 8) and taphonomy (surface condition of foraminiferal tests). Foraminiferal taxonomy constrains provenance through ecology, and taphonomy can determine residence time and transport history. For foraminifera and cave chalk analysis, we subsampled 5 cm 3 samples from each layer, wet-sieved them at 463 mm, dried at 25°C and examined them under a binocular microscope. Sieved samples were dry split to obtain counts of B300 foraminifera per sample 58 . For each sample, the total number of foraminifera present in 5 cm 3 (total concentration) was calculated, as was the percentage of each species present (abundance). Foraminiferal taxonomy followed Loeblich and Tappan 59 . Foraminifera were further categorized according to the taphonomic (surface) condition of the tests as defined by Pilarczyk and Reinhardt 51 , where pristine individuals are those that are taphonomically unaltered; abraded individuals are those that are edge rounded and corroded; and fragmented individuals are those that are broken with angular edges (Supplementary Fig. 3). In each sample, we documented the presence or absence of cave chalk. Where present, we categorized the surface condition of individual cave chalk fragments as either pristine or weathered. Pristine fragments are small, fragile and have a delicate floret structure; weathered fragments are larger and more rounded, with no floret structure.
Data availability. Data and modelling codes that have contributed from the reported results are available from the corresponding author at request.