Geological evidence of an unreported historical Chilean tsunami reveals more frequent inundation

Assessing tsunami hazards commonly relies on historical accounts of past inundations, but such chronicles may be biased by temporal gaps due to historical circumstances. As a possible example, the lack of reports of tsunami inundation from the 1737 south-central Chile earthquake has been attributed to either civil unrest or a small tsunami due to deep fault slip below land. Here we conduct sedimentological and diatom analyses of tidal marsh sediments within the 1737 rupture area and find evidence for a locally-sourced tsunami consistent in age with this event. The evidence is a laterally-extensive sand sheet coincident with abrupt, decimetric subsidence. Coupled dislocation-tsunami models place the causative fault slip mostly offshore rather than below land. Whether associated or not with the 1737 earthquake, our findings reduce the average recurrence interval of tsunami inundation derived from historical records alone, highlighting the importance of combining geological and historical records in tsunami hazard assessment. Sand sheet deposits in a tidal marsh at Chaihuín, Chile, suggest that a great earthquake in 1737 did produce a tsunami despite the lack of historical records for one. The area may be more prone to tsunami inundation than previously believed.

I nadequate anticipation of the magnitude of great earthquakes in Sumatra (2004), Chile (2010) and Japan (2011) stemmed from over-reliance on short historical records that failed to account for variability in earthquake size, rupture style, tsunamigenesis and the existence of supercycles [1][2][3] . Moreover, even where long written histories exist, individual events may be missing due to failures in reporting 4 or loss of documents in times of instability or crisis (e.g. late fourteenth/sixteenth century Japan) 5,6 , or there may be periods where only low-quality records exist containing errors 5 . Although with other limitations, geological records are free from these problems, and it is therefore imperative that we supplement historical data with geologic records, in order to obtain robust long-term patterns to inform seismic and tsunami hazard assessments.
We show that this is the case for the area affected by the 1960 Chile earthquake (magnitude 9.5), in the 1000-km-long southernmost portion of the subduction zone formed between the Nazca and South America plates (Fig. 1). Besides the 1960 earthquake and other lesser events, this area experienced three great earthquakes in the historical past (since the midsixteenth century): in 1575, 1737, and 1837. Historical and geological evidence in the form of tsunami deposits in coastal environments [7][8][9][10][11][12] , and shaking-induced turbidites in lakes 13,14 and fjords 15 , suggest that the 1575 earthquake was similar in size and extent to that in 1960, and that both were much larger than those in 1737 and 1837. The combined evidence also suggests that large tsunamis accompanied the 1575, 1837 and 1960 earthquakes but did not in 1737. To explain the complete absence of tsunami reports in 1737 and the along-strike distribution of building damage, Cisternas et al. 16 proposed a deep interplate rupture limited to the northern half of the 1960 rupture area, and restricted to beneath land, with limited capacity for tsunami generation. However, the lack of chronicles of a tsunami could also be attributed to uprisings that had driven settlers from most of the colonial outposts in the area 9,16 .
Given the limited eyewitness reports for pre-1960 south Chilean earthquakes and tsunamis, their hypothesised rupture areas and proposed mechanisms require testing with geological evidence from coastal sites. Between the long palaeoseismic records from Tirúa 7,8,17,18 (38.3°S) and Maullín 9 (41.6°S), spatially there is a gap in coastal geological evidence (Fig. 1), in which the effects of pre-1960 events are unknown. This paper addresses this gap by presenting diatom and sedimentological evidence for historical seismic events from a tidal marsh at Chaihuín, near Valdivia, close to the region of maximum 1960 slip. We aim to (1) identify and determine the timing of multiple earthquake and tsunami events from the sedimentary record; (2) use a diatom transfer function to quantify vertical coseismic deformation; and (3) test hypotheses for pre-1960 rupture areas derived from limited historical and geological records.
Tidal observations before and after the 1960 earthquake suggest 0.7 ± 0.4 m of coseismic subsidence occurred at Chaihuín (39.95°S, 73.58°W; Fig. 2) in this event 19 , and sedimentological investigations confirm subsidence 20 . With respect to pre-1960 earthquakes, to date, the closest sites to Chaihuín providing geological evidence for the 1575 earthquake are~180 km to the Fig. 1 Geographical context of the study location and evidence of historical earthquakes in south-central Chile. Location of study area a within the regional tectonic setting and b within south-central Chile relative to sites referred to in the text and relative to rupture zones of the largest historical earthquakes. Geological evidence comes from Cisternas et al. 9 15 ; historical evidence from Cisternas et al. 9,16 . Rupture lengths of historic earthquakes compiled from Cisternas et al. 16 , Dura et al. 7 and Melnick et al. 54 (dashed lines indicate where rupture extents are uncertain). Plate motions from Angermann et al. 55 . Rupture zone maps adapted from original figure drawn by Marco Cisternas. north (Tirúa 7,8,17,18 ) and south (Río Maullín 9 ), while the 1737 and 1837 events are not recorded anywhere in tidal marsh stratigraphy 9,11 , and their preservation in the geological record is limited to the south-central Chilean Lake District (39-40°S) 14 and coastal lowlands on Chiloé 10 . Chaihuín lies within the proposed rupture zone of the 1575 earthquake and tens of kilometres north of the proposed 1837 rupture area (Fig. 1). If the 1737 earthquake ruptured the northern half of the 1960 rupture area 16 , Chaihuín occupies a key location for searching for its geological signature and for delimiting the northernmost extent of the 1837 rupture. Here we present geological evidence for a previously unreported tsunami in south-central Chile, which critically reduces the Fig. 2 Location of tidal marsh and coring sites at Chaihuín. The study site is situated~20 km southwest of Valdivia (a), occupying a sheltered embayment behind a large sand spit (~2 km long,~500 m at its widest), which provides shelter from the Pacific Ocean that has promoted the development of tidal marshes and the accumulation of fine-grained sediments behind (b). The elevation distribution of the spit is extracted from the 1 m LiDAR Digital Terrain Model (d). The tidal range is low (mean higher high water [MHHW] = 0.526 m above mean sea level) and the marsh is close to the sea (<1.5 km), meaning modern diatoms should respond to tidal flooding. The sheltered location also increases the potential for preservation of evidence of seismic landlevel changes and tsunami inundation. Coring locations shown on b, c-all cores (yellow dots) are presented in Aedo et al. 22 , transects X-X' and Y-Y' are presented in Fig. 3

Results and discussion
The Chaihuín stratigraphy. Core transects (Fig. 2b) reveal three sand layers, intercalated between herbaceous peats, that are laterally extensive over 600 m across the marsh (Fig. 3a). In all cases, the sand layers have sharp lower contacts and transitional upper contacts. Ten accelerator mass spectrometric (AMS) radiocarbon dates modelled using a Bayesian phased sequence model provide the chronology (Fig. 3c and Supplementary Table 1). The age of plant macrofossils immediately beneath the upper layer, sand A, are consistent with burial by the 1960 tsunami. The age model places the deposition of the middle sand B at 1600-1820 and lower layer, sand C, at 1486-1616 CE. The calibrated age ranges for sands B and C are reasonably broad due to plateaux in the radiocarbon calibration curve, which affect dates from the seventeenth to twentieth centuries 21 .
The sedimentology and mineralogical signatures of the sand sheets are described in detail elsewhere based on over 100 handdriven cores 22 and summarised in Supplementary Discussion; here we analyse diatoms in three representative cores and present reconstructions of marsh surface elevation change over time from a diatom-based transfer function ( Fig. 4 and Supplementary Data 1). From diatom analysis of the three cores, we identified 170 species indicative of differing tolerances to tidal inundation. Only 14 species were absent from a previously published modern training set that includes 29 samples from Chaihuín 20 , and 9 of these species constituted <1% of any one sample. Fallacia ny was the only species absent from the modern samples that constituted >2% of any sample (comprising 4-5% in 2 non-sand samples).
The laterally extensive uppermost coarse to medium-grained sand sheet (A) is mid grey, varies in thickness between 1 and 19 cm, has a median grain size of 0.49 mm and is upwards fining (0.27-0.71 mm) in 61 cores (80% of those in which A is preserved, massive in the others). The marsh grades steeply into freshwater scrub, and there is no sand unit in cores just above the high marsh limit. There is an abrupt contact between the sand and dark brown silty herbaceous peat below, which contains plant material including below-ground stems (rhizomes) of Scirpus americanus. In many cores, there are rip-up clasts (~2 cm) of peat encased in the sand sheet, as well as vegetation rooted in the peat below. The peat below the sand sheet contains a diatom assemblage that is almost entirely composed of species found on the contemporary high marsh above mean higher high water (MHHW) (e.g. Eunotia praerupta, Nitzschia acidoclinata), with higher elevation optima than the diatoms found in the herbaceous peat above the sand unit (e.g. Rhopalodia constricta) (Fig. 4a). The overlying peat also contains low, albeit important, percentages (5-24%) of taxa with elevation optima below MHHW. By contrast to the peats, sand A is dominated by species with lower elevation optima (59-72% of the total assemblage have optima below MHHW), including Achnanthes reversa and Planothidium delicatulum.
The middle brown-grey to dark grey mica-rich coarse to medium-grained sand sheet (B) is similarly laterally extensive across the entire marsh, varying in thickness between 2 and 32 cm. It has a median grain size of 0.47 mm and is upwardsfining (0.38-0.68 mm) in 31 cores (50% of those in which B is preserved, massive in others), but rip-up clasts of peat were only occasionally observed. In some cases, we observe a 2-4-cm-thick cap of horizontally bedded detrital plant fragments and wood at the top of the sand layer. The sand sheet abruptly overlays a redbrown to dark brown silty herbaceous peat with variable silt content and humification. Humidophila contenta dominates the diatom assemblage in the peat below sand B (up to 37% of the assemblage) and is also present in the peat overlying the sand sheet, which remains dominated by species with elevation optima above MHHW. In the core from the lowest contemporary marsh elevation (CN14/5, Fig. 4a), there is an increase in low marsh diatom species (elevation optima below MHHW) above the sand compared to below (e.g. A. reversa, P. delicatulum). Diatom assemblages are relatively consistent across the five samples from the sand unit, with 54-76% of the assemblages being species with elevation optima below MHHW, including A. reversa, Fallacia tenera and P. delicatulum.
A third sand deposit (C) is found in 16 cores at the southern end of the marsh, although still traceable over 200 m and across most cores that penetrated deep enough to potentially sample sand C. The deposit is a dark grey fine to medium-grained massive sand (median grain size 0.25 mm, range 0.22-0.29 mm), with a maximum thickness of 51 cm and contains occasional ripup clasts from the buried organic unit below encased in the sand. The basal contact is abrupt, with the sand overlying a brown clayey silt with occasional herbaceous plant remains, humified organic matter and woody plant material. The organic horizon below sand C contains more diatom species typically found at lower elevations in the tidal frame than the peats below A and B (Fig. 4a). There is also a change in species composition approaching the top of the peat, with abundances of Opephora pacifica and Pseudostaurosira perminuta decreasing and H. contenta and E. perpusilla increasing from the base to top of the peat below sand C. Also in contrast to the other two buried organic deposits, there is a change in species composition approaching the top of the peat and samples immediately above and below sand unit C have very similar diatom assemblages, dominated by H. contenta and E. perpusilla. Diatom preservation in the sand unit was very poor, and it was not possible to obtain representative counts from this unit.
Brown silty herbaceous peats separate the three sand sheets, deposited intertidally on the basis of their diatom composition. In addition to the relative variations in freshwater and brackish diatom composition of peats described above, the peat units also vary in their degree of humification. While peats below sands A and C contain humified organic matter, the peat below sand B is unhumified. Additionally, two layers of highly humified black peat were observed immediately above and below sand A in low marsh cores from the southwest of the marsh, varying in thickness between 1 and 15 cm.
Evidence for a locally sourced tsunami. We interpret all three sand sheets as being deposited by locally sourced tsunamis, rather than far-field tsunamis or non-seismic processes (e.g. storms, river floods or aeolian processes). This is based primarily on coincident land deformation, and also upon their lateral extent, diatom composition, and sedimentological signatures. Dealing first with the latter lines of reasoning, sands A and B are not only dominated by marine sublittoral and epipsammic diatom species but also contain substantial numbers of benthic silty intertidal mudflat and freshwater taxa, which also dominate the underlying peats. This is consistent with mixed diatom assemblages in tsunami deposits worldwide and indicative of tsunamis eroding, transporting and redepositing diatoms from diverse environments as they cross coastal and inland areas [23][24][25][26] . The presence of marine and tidal flat diatoms excludes deposition of sand by river flooding 25,27 , and statistical comparison of the sedimentological and mineralogical signatures of the sands with modern depositional environments, reported by Aedo et al. 22 and summarised in Supplementary Discussion, further supports a seaward rather fluvial sediment source. We observe a maximum sedimentary contribution of 12% from upstream fluvial sources (Fig. 3b) and do not observe erosional or depositional features characteristic of fluvial flood deposits, such as a high basal mud content reflective of suspended loads during the initial stages of flooding or inverse grading as energy increases 28 .
Meteorologically driven deposition of the sands, either during storm surges or other transient sea-level fluctuation events (e.g. El Niño), is discounted as the diatoms in the overlying organic units demonstrate lasting ecological change 27,29 . While a nontsunamigenic earthquake followed closely in time by a large storm surge may impact diatom assemblages in the same way, there are several further characteristics of the three sand sheets which are consistent with a tsunami origin, even though these, in themselves, are not diagnostic. These include the lateral extent (traceable across 230 m), upwards-fining grain size of sand sheets A and B, and clasts of underlying peats observed within sands A and C and occasionally within B. The absence of extreme climatic phenomena, such as hurricanes and tropical storms, in the Chaihuín area during the historic period also minimises the possibility of finding storm deposits. However, while it is recognised that the above criteria cannot be used individually to confirm tsunami deposition, it is the combination of all sedimentological and diatom evidence that we use here in support of the most compelling evidence for tsunami deposition, which comes from the accompanying abrupt land-level change. The latter rules out deposition by tsunamis sourced in the far-field, storms or aeolian processes. Evidence for coseismic land-level change. Following established criteria 30,31 , we use the sedimentary and diatom evidence to propose that the Chaihuín sequence records three earthquake events, associated with vertical coseismic deformation and tsunami deposition. Diatom assemblages from immediately below sand layers A and B are characterised by species with higher elevation preferences than those found immediately above the sands, suggesting decreases in marsh surface elevation consistent with coseismic subsidence (Fig. 4). Diatom assemblages show minimal change across sand layer C; instead a transition occurs prior to event C whereby species with lower elevation preferences are replaced by those with higher elevation preferences, indicating net emergence prior to event C followed by minimal coseismic subsidence.
The transfer function reconstructs 0.35 ± 0.42 m of subsidence occurred in event A, which local testimony and radiocarbon dating confirm to be the 1960 earthquake. Compared to our previous estimate for this event 20 , refining the transfer function method and expanding the modern training set here, reduces the uncertainty by 0.26 m. Reconstructed subsidence agrees with observations of 0.7 ± 0.4 m 19 . By contrast, the transfer function reconstructs very minor subsidence of 0.10 ± 0.36 m occurred in event C, but this needs confirmation from analyses of additional cores.
The transfer function predicts that coseismic subsidence occurred in event B, with reconstructions varying between 0.10 ± 0.33 and 0.52 ± 0.39 m, and averaging 0.22 ± 0.38 m (Fig. 4d). While this is close to the detection limit of coseismic land-level change 30 and the error term is large compared to the amount of deformation, we interpret event B as being associated with net submergence for two reasons. First, changes in diatominferred marsh elevations between pre-and post-earthquake samples are greater than other sample-to-sample changes. Second, all nine reconstructions, regardless of core location or transfer function approach, indicate submergence rather than a mixture of submergence and emergence (Fig. 4d).
Linking the geologic and historical records. Despite the broad modelled age ranges for events B and C of 1600-1820 and 1486-1616 CE, respectively, each range only includes one historically reported earthquake. If the historical catalogue is complete, sands B and C represent tsunamis accompanying the 1737 and 1575 earthquakes, respectively. Although other great tsunamigenic earthquakes occurred in the time range of event B (1657, 1730, 1751), their rupture areas have been placed much further north 8,32 and therefore are very unlikely sources for the observed deformation. Age ranges do not include 1837; therefore, absence of evidence for this earthquake at Chaihuín supports the chronicle-based interpretation that the 1837 rupture area lies further south 11,16 . The preservation of turbidites from 1837 at sites to the north of Chaihuín 14 is consistent with observations of earthquake-triggered turbidites some distance outside the rupture zone, as observed for the M w 8.8 2010 Maule earthquake 14 .
Implications for the rupture depth in 1737. The Chaihuín record provides the first evidence for crustal deformation during the 1737 earthquake and the first evidence for the earthquake being tsunamigenic. While the nearshore bathymetry and orientation of the coastline may amplify tsunami inundation and the abundant sediment source may enhance the potential for evidence creation during even moderate tsunamis, the direction of land-level change at Chaihuín (subsidence) calls for reconsideration of the associated rupture depth. While correlation with evidence of shaking-induced turbidites from Calafquén and Riñihue lakes 14 , along with the absence of a 1737 event in sedimentary records from Río Maullín and Chucalén to the south 9,11 , supports the hypothesis that a smaller section of the plate interface ruptured in 1737 (between 39 and 41°S) than in 1960 and 1575 14 , the Chaihuín record also forms an important constraint on the depth of local slip in 1737.
By combining deformation and tsunami modelling, we show that our evidence of coastal subsidence and tsunami inundation at Chaihuín is better explained by offshore, shallow megathrust slip rather than by deeper slip below land as previously suggested 16 ( Fig. 5 and Supplementary Fig. 1). This is demonstrated by a simple numerical experiment designed to find the most likely depth range of the causative earthquake rupture that can explain the coastal subsidence inferred at Chaihuín and also the tsunami inundation.
Our numerical approach (see also "Methods") leverages the sensitivity of the deformation sign (uplift or subsidence) and tsunami size at the Chaihuín coast to the depth of megathrust slip 33 (Fig. 5). An earthquake rupture with maximum slip at 33 km fault depth (Fig. 5a, red model), as previously inferred from historical records 16 , will result in coastal uplift and a relatively small tsunami. Instead, if the rupture occurs offshore (Fig. 5a, blue model), the deformation will result in coastal subsidence and a much larger tsunami. From a systematic analysis in which the hypothetical rupture models are shifted horizontally in the trench-normal direction or vertically in the depth direction (Fig. 5a, upper panel), we conclude that subsidence at the Chaihuín coast could only be produced by ruptures placed mainly offshore, at average megathrust depths shallower than 20 km (Fig. 5b, downward triangles). Deeper ruptures will produce coastal uplift and consequent smaller tsunamis (Fig. 5b). The same conclusion is reached by varying the rupture width with fixed updip and downdip limits (Supplementary Fig. 1).
Our conclusions are independent of the use of a normalised unit displacement in all models (i.e. 1 m at the centre of its corresponding bell-shaped rupture) because the opposing effects of deep versus shallow ruptures at Chaihuín are insensitive to the magnitude of slip involved and depend on its locus. The amount of slip determines the magnitude of deformation but not its sign due to the elastic response of the crust during earthquakes 34 . However, with evidence at only one location we only feel confident to constrain the depth range but not the magnitude nor along-strike extent of the causative slip. Therefore, from our numerical experiment we conclude that to produce subsidence at the Chaihuín coast, an offshore rupture likely shallower than 20 km is required as a deeper source would result in coastal uplift. This is also consistent with the inferred tsunami heights (Fig. 5b), which are larger for a shallower rupture and therefore more likely to produce inundation on land independent of the local topography. This geologically-based inference of an offshore rupture (blue curve in Fig. 5b) contrasts with the deeper rupture below land (red curve in Fig. 5b) previously inferred from historical observations alone 16 .
Implications for tsunami recurrence intervals. The average interval between the three events preserved at Chaihuín, 193 years, is shorter than the interval proposed for full segment 1960style ruptures of 270-280 years 9,11,14 . This supports the notion that the Chilean subduction zone displays a variable rupture mode, in which the size, depth, tsunamigenic potential and recurrence interval vary between earthquakes 10 . Of greatest importance, however, is the shorter average recurrence interval of tsunami inundation than previously reported. With the addition of the 1737 tsunami alongside previously known events in 1960, 1837 and 1575, the historical recurrence interval for tsunamis generated anywhere along the Valdivia segment of the Chilean subduction zone is reduced to 130 years. This holds even if the inferred tsunami inundation is not associated with the 1737 earthquake, but with another earthquake of similar age missed in the historical catalogue.

Conclusions
Locally sourced tsunami deposition is the favoured mechanism for the emplacement of three sand sheets at Chaihuín tidal marsh, south-central Chile. We infer this mechanism based on the lateral extent of the deposits, sedimentological composition compared to modern environments, characteristic sedimentary structures including rip-up clasts, marine diatom assemblages that exclude a fluvial source, and most importantly, coincident abrupt marsh submergence that excludes storm surge deposition.
Tidal marsh stratigraphy and changing diatom assemblages record three episodes of abrupt land-level change accompanying the tsunami deposits. Modelled ages for the youngest and oldest deposits correspond with known tsunamigenic earthquakes in 1960 and 1575 CE. The intervening event is consistent in age with In this suite of models, the rupture width and peak slip are fixed at 100 km and 1 m, respectively, and the rupture location is systematically shifted horizontally in the trench-normal direction to represent ruptures at different depths. b Summary plot showing the modelled coastal uplift (left vertical axis) and tsunami runup (right vertical axis) predicted by the suite of models. Note that coastal subsidence can only be produced by offshore ruptures, with slip shallower than~20 km. Ruptures deeper than this produce uplift at the coast. This opposing pattern of coastal deformation between shallow versus deeper ruptures is insensitive to how much slip is prescribed at the fault. Supplementary Fig. 1 shows the results for two different suite of models, in which the rupture width varies by fixing the updip (Supplementary Fig. 1a) and downdip ( Supplementary Fig. 1b) limits.
the 1737 earthquake, for which no tsunami has previously been reported. The lack of evidence for land-level change or tsunami inundation consistent with a historically documented earthquake and tsunami in 1837 confirms Chaihuín lies to the north of the 1837 rupture area.
Our identification of the first evidence for a tsunami accompanying the 1737 earthquake reveals tsunami inundation has occurred more frequently than previously thought along the Valdivia segment of the Chilean subduction zone. Furthermore, by using coupled deformation-tsunami modelling, guided by the new geological evidence, we show that the 1737 earthquake ruptured mainly offshore, at fault depths much shallower than previously proposed from historical records alone. These results provide further support for the spatiotemporal variability in ruptures within a supercycle between the 1575 and 1960 M w > 9 earthquakes.
This study highlights that historical records may not provide complete documentation of the occurrence and characteristics of earthquakes and tsunamis occurring during the historical period. Geological evidence is essential for not only extending earthquake and tsunami histories into the past, but for verifying and supplementing historical records. This has global significance anywhere where there could be issues of absent or incomplete reporting of events, or loss of records, particularly in times of local or national crisis.

Methods
Tidal marsh sediments are excellent recorders of relative sea-level changes 35 , and in regions of palaeoseismicity, sequences of organic soils interbedded with minerogenic units often reflect vertical land-level changes that occur both coseismically during megathrust earthquakes and through intervening interseismic periods 36,37 . Diatoms, siliceous microfossils incorporated within these tidal marsh sediments, are used to quantify vertical land-level changes associated with great earthquakes due to the close control of salinity on their distribution in intertidal environments, as well as their high preservation potential in coastal sediments 38 . Due to their robust valves, diatoms are also widely used to determine the provenance of tsunami sediments and changing flow conditions during tsunamis 23,[38][39][40] . In this way, diatoms have been used to reconstruct the history of earthquakes and tsunami worldwide over multi-millennial timescales 25,31,[36][37][38] , providing long-term perspectives on the recurrence, magnitude and variability of earthquakes occurring along subduction zones.
Lithostratigraphy. We used marsh front exposures and 130 hand-driven gouge and Russian cores (collected between 2013 and 2018) to reveal the stratigraphy. We sought evidence for tsunami inundation in the stratigraphic record, but as sand sheets can be deposited by both seismic and a range of non-seismic sedimentary, fluvial, oceanographic or atmospheric processes, particularly storm surges 29,31 , we employed transects of cores (Fig. 2b) to assess the lateral extent and continuity of sand sheets. We interpret sand sheets as being of seismic origin where these abruptly overlay organic peats in cores correlated over 10s to 100s of metres, with diatom, chronological and sedimentological analyses providing further lines of evidence 30,31 .
Sedimentological and mineralogical analyses. Methods used to characterise the sediments according to grain size and mineralogy are detailed in Aedo et al. 22 . To summarise briefly here, 25 subsamples (50-140 cm 3 ) were taken from the sand sheets recovered from cores to analyse the granulometry and mineralogy. Granulometric analysis followed the particle settling velocity method 41 , and central tendency and mean dispersion were calculated using GRANPLOT 42 . The mineral composition was estimated based on a modal count with a binocular microscope of 200 grains. In order to establish the origin of the sands, core sand samples were statistically compared to 22 samples taken from a series of modern environments (beach, sand bar, dunes, estuary and the outer bay, Fig. 2b). Discriminant analysis was performed using the modern samples and four granulometric parameters (mean, selection, skewness, normalised kurtosis) as predictor variables, using the Statgraphics Centurion XVI software version 16.2.04. Finally, we also scanned selected core sections using a Geotek X-ray core imaging system at Durham University. The scanner operated at 125 kV and we used Fiji 43 to visualise the data.

Diatom analysis.
A key criteria in assessment of tsunami versus non-seismic sources of sand sheets is coincident vertical land-level change 30,31 . We therefore analysed fossil cores for diatoms due to their well-documented utility in relative sea-level reconstruction, which is based on the relationship between species and elevation 35 . Due to varying tolerances to salinity, the frequency and duration of tidal inundation and sediment type, diatom assemblages show distinct vertical zonation across a tidal marsh. We discuss these contemporary relationships between species and marsh surface elevation in Hocking et al. 20 and use the modelled contemporary relationships to quantitatively relate the changes seen in fossil diatom assemblages to past changes in the elevation of the marsh surface via a transfer function. We prepared samples for diatom analysis following standard methods 40 and counted a minimum of 250 diatom valves per sample.
In the transfer function, we use a regional south-central Chile modern diatom training set (Supplementary Data 1), updated from Hocking et al. 20 . The training set consists of 148 modern samples from 6 marshes in south-central Chile (Queule, Río Lingue, Isla del Rey/Valdivia, Chaihuín, Quilo and Estero Guillingo; Fig. 1). We use a regional modern training set due to the superiority of the regional southcentral Chile model over sub-regional and local models, both in terms of agreement with observations in 1960 and in providing fossil samples with the closest modern analogues 20 . This also follows studies elsewhere that favour regional models, particularly as the difference in age between modern and fossil samples increases and the present day marsh may no longer capture the full environmental conditions potentially present in older records [44][45][46] . While modern sample sites span~280 km, the training set does include 29 local samples from Chaihuín, as previous studies have also highlighted the importance of including samples from the local area where possible 20,44 . Since publication of the modern training set in 2017, we have collected further tidal data to improve the tidal modelling and calculation of modern sample elevations, as well as updated nomenclature. These changes have resulted in a minor change in transfer function performance and required that modern tidal flat samples and those from elevations above a standardised water level index (SWLI) of 350 be excluded from the training set (see Hocking et al. 20 for explanation of standardisation of sample elevations to account for tidal range differences between sample sites). Outside of this threshold range, the transfer function does not do well at predicting elevation. The updated modern training set contains 7 samples that were not included in 2017 and excludes a further 33 samples, which were previously included.
Development of the transfer function and evaluation of transfer function performance builds upon Hocking et al. 20 . We follow the same approach of using weighted-averaging partial least squares (WA-PLS) regression (due to the unimodal nature of the diatom data) with bootstrapping cross-validation (1000 cycles), and assessed performance by using the Modern Analogue Technique to calculate the similarity between modern and fossil assemblages (reported as minimum dissimilarity coefficients). We also run a locally weighted-weighted averaging (LW-WA) transfer function model, with cross-validation and both classical and inverse deshrinking, in order to compare reconstructions derived by different transfer function methods. The three models used in the transfer function The output of the transfer function is the palaeomarsh surface elevation (PMSE) associated with each fossil sample. To quantify coseismic deformation, we subtract the reconstructed pre-earthquake PMSE from the post-earthquake PMSE and correct for sediment accumulation, including tsunami deposition. The associated 95.4% (2σ) uncertainty term is the square root of the sum of the squared samplespecific standard errors for pre-and post-earthquake samples. In order to give an overall estimate of deformation for each earthquake, we average transfer function outputs from all cores and model approaches, rather than relying on a single method.
Chronology. We assigned ages to sand layers by using AMS radiocarbon dating of herbaceous plant macrofossils found immediately above and below the sands. We selected above ground parts of terrestrial plants that were horizontally bedded within the sediment matrix wherever possible. Samples were analysed at high precision using multiple graphite targets, due to the potential problems for radiocarbon dating and calibrating ages from the past 400 years associated with plateaux in the radiocarbon curve 21,47 . We calibrated dates using the SHCal20 calibration curve 48 or the post-bomb atmospheric southern hemisphere calibration curve 49 for the samples exceeding 100% modern carbon. We report 2σ age ranges in years CE. We used a Bayesian phased sequence model within OxCal (version 4.4) to refine the calibration solutions; this constrains the relative order and grouping of events 50 . The model includes 10 macrofossil ages; one charcoal sample was excluded due to the anomalous old age, which was out of stratigraphic order (Supplementary Table 1).
Rupture and tsunami modelling. We modelled coastal deformation at Chaihuín due to megathrust slip using a three-dimensional (3D) dislocation model in a uniform elastic half-space with an assumed Poisson's ratio of 0.25. We used an extensively benchmarked computer code 51 that numerically integrates the pointsource dislocation solution (Green's function) 34 over a 3D megathrust and yields displacement values on the surface of the model. Given that we are interested in constraining the location of the rupture in the dip direction rather than in the strike direction, we constructed a megathrust that is uniform along strike and long enough to avoid 3D or edge effects. The megathrust geometry was constructed by discretizing slab data 52 at the latitude of Chaihuín into small triangular fault elements, each representing individual point-sources. Deformation at the surface was computed by the sum of the contributions from all the triangular sources. All tested earthquake rupture models consider a bell-shaped slip distribution in the dip direction with a peak slip of 1 m at its centre, which tapers updip and downdip 51 . The peak slip is fixed at unity because the sign of coastal deformation (either subsidence or uplift) depends only on the rupture depth. Constraining the slip magnitude or how it varied along-strike would require additional evidence at other sites along the south-central Chile coast. All models assume that the direction of coseismic slip (rake) is roughly 90°to represent pure dip-slip faulting.
All tsunami simulations were computed by using a finite-difference method on the actual bathymetry offshore Chaihuín. We used the well-validated numerical model COMCOT, which solves the linear shallow water equation (LSWE) and nonlinear shallow water equation using a leap-frog scheme on a staggered and nested grid system 53 . We assumed the initial sea-surface elevation to be equal to the vertical seafloor deformation due to earthquake faulting of each rupture model. The two-dimensional bathymetric profile, which we uniformly extended along strike, was constructed from global gridded bathymetric data (www.gebco.net), from which a two-level nested grid system was used. The first grid level with 500 m spatial resolution was used in the source region and in deep water, where the LSWE were considered. In the coastal area of Chaihuín, we upsampled the GEBCO data by a factor of 5 to obtain a finer grid of~100 m resolution for shallow water propagation and runup computation. Even though using a realistic, highresolution nearshore bathymetry will change the details of the results (i.e. modelled tsunami runup values in Fig. 5b), our simplified coarse bathymetric model is sufficient for the purpose of the numerical experiment described above. We assumed an open radiation boundary condition to the outer sea and a moving boundary on the coast. The computation time step for all grids was set to satisfy the Courant-Friedrichs-Lewy stability condition of the finite difference method. All tsunami simulations were run for 1 h, which we checked was long enough to obtain the largest runup at the coast (in this case associated with the leading wave).

Data availability
Diatom, sedimentology, stratigraphy and radiocarbon data sets are available on the Figshare repository with the identifier: https://doi.org/10.6084/m9.figshare.16617241. Data sets used for tsunami modelling are freely available online or from linked references in the "Methods" section (COMCOT tsunami code and dislocation model 53 ; SLAB geometry 52 ; GEBCO bathymetry from www.gebco.net). LiDAR data were donated by Forestal Arauco to the CYCLO project under a confidentiality agreement and may be obtained from the corresponding author upon reasonable request.

Code availability
The COMCOT tsunami code is available from Wang 53 and the dislocation modelling code as described in Gao et al. 51 was provided by Kelin Wang. No new computer code was written in preparing the paper.