Modern diatom assemblages from Chilean tidal marshes and their application for quantifying deformation during past great earthquakes

Tidal marsh sediments from south‐central Chile provide evidence for multiple great earthquakes. Diatom transfer functions, statistical models of the relationship between species preserved in the sediment and elevation, provide quantitative estimates of coseismic vertical land‐level change associated with individual earthquakes. However, in south‐central Chile, our ability to quantify land‐level change is currently limited by a lack of understanding of the environmental variables controlling the distribution of diatoms, an essential prerequisite for converting variations in fossil diatom assemblages into quantitative estimates of past elevation changes. We present a new modern diatom dataset for the region and explore the implications of the scale of the dataset used in transfer function models on the reconstructions of land‐level change. Modern training sets containing samples from a regional scale are superior to sub‐regional and local‐scale training sets, providing closer estimates for known deformation during the great 1960 Chilean earthquake, a higher proportion of good modern analogues and uncertainty terms up to 42% smaller than previously published reconstructions. Copyright © 2017 The Authors. Journal of Quaternary Science Published by John Wiley & Sons Ltd.


Introduction
Great earthquakes (magnitude >8) occur along subduction zones, with cycles of strain accumulation and release resulting in gradual interseismic deformation punctuated by episodes of abrupt coseismic deformation. The litho-and biostratigraphy of tidal marshes may preserve a record of the vertical component of deformation (e.g. Atwater, 1987;Zong et al., 2003), and have been used to quantify the magnitude of coseismic uplift or subsidence and reconstruct the occurrence of earthquakes over multi-millennial timescales (e.g. Combellick, 1991;Nelson et al., 1996;Shennan et al., 1996;Sawai et al., 2004;Hamilton and Shennan, 2005). Such palaeoseismic records provide long-term perspectives on the recurrence, magnitude and variability of earthquakes occurring along subduction zones and contribute to seismic hazard assessments (e.g. Wesson et al., 2007;Mueller et al., 2015). Sedimentary and biostratigraphic evidence attests to the occurrence of past episodes of coseismic deformation and tsunami inundation along the Chilean subduction zone (Atwater et al., 1992;Bartsch-Winkler and Schmoll, 1993;Cisternas et al., 2005;Nelson et al., 2009;Garrett et al., 2013Garrett et al., , 2015Dura et al., 2015). However, incomplete understanding of present-day tidal marsh environments currently limits our ability to accurately quantify land-level changes. Quantitative methods for reconstructing land-level change are based on uniformitarianism and the assumption that the relationships between contemporary organisms and their environment can be used to reconstruct past environments (Birks, 1995;Birks et al., 2010). Diatoms (unicellular algae) are widely used for this purpose due to their vertical zonation in intertidal environments, as well as their high preservation potential in coastal sediments (for full reviews see Pilarczyk et al., 2014 andDura et al., 2016). In Chile, modern diatom data are available only from a limited number of sites (Nelson et al., 2009;Garrett et al., 2013;Sawai et al., 2017). Consequently, as we attempt to reconstruct palaeoseismic deformation at new sites or further back in time, our models may be limited by not having a modern-day equivalent of the past environment we are trying to reconstruct. Particularly in environments which have undergone rapid decimetre-to metre-scale land-level change during an earthquake, the present-day marsh may not provide a good modern analogue for an entire fossil sequence. In addition, reconstructions of past land-level change can also vary depending on the spatial scale of investigation, and whether fossil diatom assemblages are compared to modern samples collected from the same marsh, or from a regional ensemble of samples (Horton and Edwards, 2005;Watcham et al., 2013).
This paper seeks to improve understanding of modern intertidal diatom distributions in south-central Chile, the region that produced the world's largest recorded earthquake, the magnitude 9.5 Valdivia earthquake in 1960. This will allow more accurate and precise quantification of palaeoseismic deformation using a transfer function approach. We aim to summarize modern diatom distributions from four new tidal marsh sites, combine these data with previously published diatom distributions from two sites on Isla Chilo e (Garrett et al., 2013) and investigate the relationship between diatom assemblages and elevation. Following the approach applied by Watcham et al. (2013) in Alaska, we create a series of modern datasets structured by spatial scale, and investigate the implications of using modern data from local, sub-regional and regional scales on model performance and coseismic deformation reconstructions. Finally, we assess the viability of using our expanded modern training set for reconstructing deformation associated with past great earthquakes.

Tidal marsh sites
Tidal marsh sediments are excellent recorders of changes in relative sea level ; in a region of palaeoseismicity, sequences of organic soils interbedded with minerogenic units reflect vertical deformation throughout the earthquake deformation cycle. The sedimentary signatures of seismic land-level change are described in detail elsewhere (e.g. Atwater, 1987;Darienzo and Peterson, 1990;Shennan et al., 1996;Cisternas et al., 2005;Hamilton and Shennan, 2005;Garrett et al., 2015), along with criteria used to differentiate organic-minerogenic couplets of seismic origin from those potentially deposited by other non-seismic sedimentary, fluvial, oceanographic or atmospheric processes, particularly storms (Nelson et al., 1996;Witter et al., 2001;Shennan et al., 2016). Tidal marshes may also preserve evidence for tsunami inundation through the emplacement of laterally extensive sand sheets; the identification of these deposits provides an additional criterion for inferring a seismic origin for sedimentary couplets (Nelson et al., 1996;Witter et al., 2001;Shennan et al., 2016).
We focus our investigation on tidal marsh sites within the 1960 earthquake rupture area, along the Chilean subduction zone (37.5-46˚S; Fig. 1). The 1960 rupture caused uplift in a predominantly offshore belt and up to 2 m of subsidence on the coast (Plafker and Savage, 1970;Fig. 1). Studied south-central Chile tidal marshes were chosen for two reasons: (i) previously published works demonstrate their potential as good sedimentary recorders, especially compared to those of central Chile (May et al., 2013), and (ii) detailed observations of land-level change in 1960 are available with which we can compare our diatom-based reconstructions (Table 1).
We combine new modern diatom data from four tidal marsh sites in south-central Chile À Queule, Río Lingue, Isla del Rey and Chaihuín À with those already published from two sites on Isla Chilo e (Quilo and Estero Guillingo; Supplementary Information, Fig. S1; Garrett et al., 2013). We then use this new modern dataset to quantitatively reconstruct land-level change associated with the 1960 earthquake at five sites À Río Lingue, Chaihuín, Pucatrihue, Llico and Río Maullín (Fig. 1). We also re-quantify landlevel change at three previously investigated sites (Río Andali en, Río Tir ua and Chucal en; Garrett et al., 2013Garrett et al., , 2015 using the newly expanded modern diatom training set. All the tidal marshes sampled occupy sheltered embayments alongside estuaries, and in the case of Chaihuín, Pucatrihue and Llico, large sand spits provide further shelter from the ocean (Fig. 1B). Their sheltered locations promote the development of tidal marshes and the accumulation of fine-grained sediments, increasing the potential to preserve evidence of seismic land-level changes. All the studied marshes are connected to and are close to the sea (<2 km), meaning modern diatoms should respond to tidal flooding. Finally, none of the modern sample sites were affected by coseismic deformation or ongoing post-seismic deformation associated with the 2010 Maule earthquake further north.  Plafker and Savage (1970); precise values of deformation for individual sites are provided in

Field and laboratory methods
To obtain a training set of modern diatom assemblages, in January 2013 we sampled at low tide the top $1 cm of sediment (after removing the surface film) at regular vertical intervals (every 5 cm, with the exception of the highest marsh at Chaihuín, which we sampled at 10 cm vertical intervals). We took modern samples along transects perpendicular to the direction of tidal inundation, covering the full environmental gradient from unvegetated mudflat to freshwater marsh. Additional details on modern sampling locations and dominant marsh vegetation are included in Fig. S1. The number of samples taken per transect depended on the elevation gradient of the marsh. Marsh front exposures, pits, gouge and Russian cores provided stratigraphic evidence for past earthquakes. At all sites we employed transects of cores to assess the lateral continuity and abruptness of sedimentary couplets and the presence of tsunami-deposited sand sheets. We interpreted couplets as indicative of coseismic subsidence only where organic peat is abruptly overlain by sediment indicative of a lower elevation in cores correlated over tens to hundreds of metres (Table S1), with diatom and chronological analyses (discussed below) providing further lines of evidence (cf. Nelson et al., 1996;Shennan et al., 2016). We interpreted sand sheets as being indicative of tsunami deposition on the basis of key characteristics (e.g. abrupt lower and upper contacts, often with flattened vegetation rooted in the underlying unit, laterally extensive, landward thinning, upward fining, rip-up clasts, mixed diatom assemblages, following Morton et al., 2007;Peters and Jaffe, 2010;Engel and Br€ uckner, 2011). We levelled each modern sample and fossil core location to a benchmark and related this to mean sea level (MSL) by comparison of tidal measurements from an ultrasound tide gauge (Wesson et al., 2014) deployed at each field site with the nearest permanent tide gauge station (IOC Sea Level Station Monitoring Facility; Table 2). To allow for differences in tidal range between sites, we convert sample elevations to a standardized water level index (SWLI), whereby an SWLI value of 100 represents MSL and 200 represents mean higher high water (MHHW) (following Hamilton and Shennan, 2005): where: SWLI n is the standardized water level index of sample n, h n is the elevation of sample n (in metres), h MSL is mean sea level at the site (in metres) and h MHHW is mean higher high water at the site (in metres). We analysed modern and fossil samples for diatoms, due to their well-documented utility in relative sea-level reconstruction, which is based on the relationship between species and elevation (Nelson and Kashima, 1993;Horton, 1998, 1999;Sherrod, 1999). Due to varying tolerances to salinity, the frequency and duration of tidal inundation, and sediment type, we expect diatom assemblages to show distinct vertical zonation across a tidal marsh (Fig. S1). We Table 1. Observations of vertical land surface deformation caused by the 1960 earthquake from Plafker and Savage (1970 Plafker and Savage (1970) obtained in 1968 from (i) tides: measured difference between the preand post-earthquake extreme tides as indicated by local residents, (ii) vegetation: measured difference between the pre-and post-earthquake lower growth limit of terrestrial vegetation, and (iii) local residents: estimated change in the position of shoreline markers relative to tide levels as reported by local residents. ‡The error terms reflect the error of the method of observation and for sites marked by Ã the negative error term also includes possible local sediment compaction of 0.5 m. Starred localities are on unconsolidated deposits where some surficial subsidence may have occurred. Errors associated with the method of observation as cited by Plafker and Savage (1970): measured changes from vegetation growth limits are considered minimum estimates as pre-earthquake vegetation limits were commonly eroded shoreward. Accuracy of measurements was considered good (AE0.2 m) where both pre-and post-earthquake vegetation lines were sharply defined, and fair (AE0.4 m) where only one limit was clearly defined. Similarly, the accuracy of tidal measurements was considered good where there was reasonable certainty of pre-earthquake extreme high tides relative to fixed markers in bedrock areas, and fair where pre-earthquake tide levels could only be approximately located (Plafker and Savage, 1970). can use the contemporary relationships between species and marsh surface elevation to relate the changes we see in fossil diatom assemblages to changes in the elevation of the marsh surface in the past using a transfer function approach (see 'Transfer functions'). We prepared samples for diatom analysis following standard methods (Palmer and Abbott, 1986) and counted a minimum of 250 diatom valves per sample. We confirmed that sedimentary organic-minerogenic couplets relate to the 1960 earthquake by a combination of accelerator mass spectrometry radiocarbon dating, local testimony and elevated caesium-137 concentrations in organic soils immediately below sand layers, attributed to nuclear weapons testing that began in the early 1950s and reached its maximum radioactive fallout in 1963 (Pennington et al., 1973;Arnaud et al., 2006). For the last of these, we compared caesium-137 concentrations in the sample immediately underlying the sand layer with background concentrations in an additional deeper sample. For radiocarbon dating, we selected above-ground parts of terrestrial plants that were horizontally bedded within the sediment matrix, and calibrated dates in OxCal (version 4.2) using the postbomb atmospheric southern hemisphere calibration curve (Hua and Barbetti, 2004).
Where cores had been previously collected, described and dated (Río Andali en, Río Tir ua, Río Maullín and Chucal en), diatom assemblages were recalibrated using the new modern diatom training set.

Development
To reconstruct changes in marsh surface elevation across observed sedimentary transitions in our fossil sequences, we use a diatom-based transfer function. This is a two-stage process involving regression to establish the present-day relationship between diatom species and elevation, and then calibration to use the modelled contemporary relationship to quantify the palaeomarsh surface elevation (PMSE) associated with each fossil sample.
Detrended canonical correspondence analysis confirms there is a unimodal species response to elevation in the modern environment at all sites (environmental gradient lengths >2 standard deviations; Table 3). Therefore, we use the weighted averaging partial least squares (WA-PLS) regression technique that is recommended where speciesenvironment relationships are non-linear, with bootstrapping cross-validation to obtain sample-specific errors (ter Braak and Juggins, 1993;Birks, 1995). We select the model with the strongest and most linear correlation between observed and predicted elevation values, which at the same time minimizes the root mean square error of prediction (RMSEP) and maximum bias. As increasing the number of model components can lead to over-fitting, where the differences in model performance are small and ambiguous, we follow the principle of parsimony and choose the 'minimal adequate model', i.e. the one with fewest parameters (Crawley, 2005;Juggins and Birks, 2012). As well as considering the model performance measures, we also consider the effects on species optima as increasing components are added (cf. Wright et al., 2011). We then apply our chosen transfer function model to the fossil samples and convert the resulting PMSE estimates from SWLI units to metres by rearranging Eqn (1). 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% (2s) error term is the square root of the sum of the squared sample-specific standard errors for pre-and postearthquake samples. It is not known whether tsunami sand emplacement is associated with erosion of the pre-earthquake marsh. Therefore, reconstructions of coseismic subsidence and uplift should be regarded as minima and maxima, respectively. However, as reconstructed PMSEs of preearthquake samples used for reconstructing coseismic landlevel change and those samples immediately below are consistent at each site to within a maximum of 9 cm (2 cm at five out of eight sites), this suggests minimal pre-earthquake PMSE variation and therefore the impact on reconstructions of any potential erosion of the pre-earthquake marsh is likely to be minimal.

Evaluation
We assess the reliability of our reconstructions through (i) numerical evaluation, (ii) comparison with field observations where available and (iii) investigation of the elevation signatures of principal diatom taxa in the modern environment. We use the modern analogue technique (MAT) to quantify the similarity between each fossil sample and the modern training set (Birks, 1995). We aim to select modern training sets that provide the most fossil samples with good, or at least close, modern analogues, whereby we define 'good' and 'close' as the 5th and 20th percentiles, respectively, of dissimilarity values for the modern samples (Simpson, 2007;Watcham et al., 2013). A minimum dissimilarity coefficient (minDC) of zero indicates the fossil and its closest modern samples are exactly similar, and the coefficient increases with increasing dissimilarity (Jackson and Williams, 2004). We developed the transfer function and performed the MAT using the C2 software package (version 1.7.3, Juggins, 2011).

Scale and sample selection
There is ongoing debate about the scale of modern training sets to use in transfer function reconstruction, between local models that generally give lower sample-specific errors versus regional models that generally do a better job of providing fossil samples with modern analogues as the difference in age between modern and fossil samples increases (Gehrels et al., 2001;Horton and Edwards, 2005;Nelson et al., 2008;Woodroffe and Long, 2010;Wilson and Lamb, 2012;Watcham et al., 2013). Some studies advocate using local training sets, on the basis that using regional modern datasets reduces precision and increases variability in diatom assemblages not related to elevation (Woodroffe and Long, 2010), while others favour a regional approach, especially when reconstructing relative sea-level change over longer periods up to thousands of years (Gehrels et al., 2001;Horton and Edwards, 2005;Nelson et al., 2008;Watcham et al., 2013;Sawai et al., 2016). Here we assess the (dis)similarity of modern samples from different sites using the unconstrained ordination technique detrended correspondence analysis (DCA) and cluster analysis (performed in PRIMER; Clarke and Gorley, 2015), and then aim to add to this local versus regional debate by comparing reconstructions based on modern training sets from both scales.
Following Birks (1998), we include all species in the models. While the distribution of some species in the modern samples may not be primarily controlled by elevation, there may also be variation in fossil assemblages explained by factors other than elevation. We do, however, investigate the effect of excluding samples based on elevation.
There is a noticeable decline in the abundance of salt-tolerant species and increase in freshwater species above MHHW (200 SWLI). The species that dominate below MHHW have broad elevation tolerances reflecting either cosmopolitan distributions and salinity tolerances or reworking and redistribution of diatom valves. We acknowledge that these species will be less useful in transfer function reconstructions but do not exclude them from modern training sets as this mixture of in situ and reworked valves will also occur in fossil sediments. Overall, we see more variable and species-rich assemblages above MHHW and relatively low species diversity and the dominance of key taxa at lower elevations below MSL. Species distributions are generally similar across the modern training sets from different locations, although there are some notable exceptions. In particular, there are a number of species unique to, or found in higher abundances at, Quilo and Estero Guillingo on Isla Chilo e (e.g. Cocconeis scutellum, Diadesmis contenta, Eunotia paludosa, Fragilaria construens var. venter, Navicula atomus, Pinnularia intermedia, Sellaphora pupula). This could reflect altered nutrient availability, as the tidal flats are used for algal harvesting, or could relate to differences in sediment type or salinity. DCA confirms the similarity between samples from sites within the Araucanía (Queule) and Los Ríos regions (Río Lingue, Isla del Rey, Chaihuín), and between samples from Isla Chilo e (Fig. 3). However, clusters from Isla Chilo e and Araucanía/Los Ríos occupy discrete areas of the DCA sample plot, suggesting there are regional differences between the diatom assemblages. Therefore, we base the groupings of modern training sets on regions. From the six local training sets, we combine Queule, Río Lingue, Isla del Rey and Chaihuín to form a  (2017) sub-regional Los Ríos training set, and combine samples from Quilo and Estero Guillingo to form a sub-regional Chilo e training set. We then combine all samples to form a regional modern training set for south-central Chile. The limited spatial variability in the elevation optima of the principal diatom taxa further justifies the compilation of a combined regional modern training set (Fig. 4B,C). The elevation optima of all principal diatom species (>15% of at least one sample in the regional training set) vary by less than 25 SWLI units (one-quarter of the difference between MHHW and MSL) between sub-regional and regional training sets, with three notable exceptions of Navicula cincta, Navicula phyllepta and Cymbella pusilla with variable optima of $50 SWLI (Fig. 4C). We discuss the impact on reconstructions of these species under 'Assessing the optimum spatial scale of the modern training' and 'Do Chilean tidal marshes underestimate coseismic deformation?'.

Transfer function development
Summary statistics of WA-PLS regression show that for most of the training sets, the two-component model performs best (Table 3). For the Chilo e and regional datasets, the addition of a third component improves the RMSEP by 5%, and this is frequently used as the criterion to justify selection of a more complex three-component model (Birks, 1998). However, here we also investigate the effect of model component selection on species optima and ultimately on reconstructions (following Wright et al., 2011). Figure 4A shows that many of the elevation optima (beta coefficients) of principal diatom species (those comprising >15%) are strongly distorted when using component 3. At sites where fossil assemblages are dominated by species such as Diadesmis contenta, whose optima are significantly altered by the addition of a third component (e.g. Pucatrihue and Río Maullín), we see that reconstructions of land-level change are consistent across different spatial scales using a two-component model, but differ when we apply the three-component model (Table 5; Fig. 4c). This latter variation is therefore likely to be a function of assemblage-specific optima updates distorting reconstructions using the third component, and we therefore cannot justify the application of component 3 solely on the grounds of improved model performance statistics.
Likewise, we also justify using two-component over onecomponent models on the basis that at some sites (e.g. Río Lingue, Chaihuín, Chucal en) there appears to be real variation, rather than the apparent variation described above, between fossil reconstructions using training sets from different scales (Table 5). Variations in the reconstructions are therefore likely to be a result of improved modern analogues rather than a shifting of species optima. We are therefore justified in using the two-component model, which has the advantage over the one-component model of improving model performance (reducing error and biases in predicted values), and being less susceptible to the 'edge effects' of weighted averaging that can lead to non-linear distortions at the end of particularly long environmental gradients (ter Braak and Juggins, 1993;Juggins and Birks, 2012). We do, however, still recommend viewing the reconstructions with caution where the fossil records are dominated by cosmopolitan diatom species with an inconsistent elevation signature across different modern tidal marshes, e.g. Navicula cincta and Navicula phyllepta (Fig. 4B).
We chose to exclude samples below SWLI values of 125, as below this threshold the transfer function does not do well at predicting elevation (Fig. 5), and the implications of this are discussed below ('Do Chilean tidal marshes underestimate coseismic deformation?'). Above this threshold, there is a linear correlation between elevations predicted by the second component of the regional transfer function model and observations, and the even spread of residuals above and below zero confirms predictions do not suffer from bias that can be associated with WA-PLS regression ( Fig. 5; Birks, 1998). The pruned modern training sets consist of between 29 and 176 samples, exceeding the minimum of 25 that is required to give sample-specific error terms on PMSE reconstructions (Telford and Birks, 2011).

Chronology
Before running the selected transfer function models, we first confirm the sedimentary couplets in the fossil record relate to the 1960 earthquake. Radiocarbon ages for Río Lingue, Chaihuín, Pucatrihue and Llico are reported in Table 4; for Río Maullín see Cisternas et al. (2005) and for Chucal en see Garrett et al. (2015). Above-background caesium-137 concentrations in the buried marsh surface  confirm that sand layers at Río Andali en and Río Tir ua relate to the tsunami associated with the 1960 earthquake (Garrett et al., 2013).
Assessing the optimum spatial scale of the modern training set We investigate how reconstructions differ when we use modern training sets from varying spatial scales, as well as how this affects model performance (Table 5). We found differences between models at each site of between 0.05 and 0.5 m in reconstructed mean marsh surface elevation change associated with the 1960 earthquake (using WA-PLS component 2). At some sites (e.g. Río Tir ua, Río Lingue, Pucatrihue and Río Maullín), the variation in reconstructions between models equates to <30% of the maximum coseismic deformation reconstructed for that site. By contrast, at Llico there is a difference of 0.43 m between the mean reconstruction from the sub-regional Los Ríos and regional scale models (70% of maximum coseismic land-level change for this site). This is probably due to the dominance of two species in the fossil record whose elevation optima are sensitive to training set selection (the elevation optima of Luticola mutica and Pinnularia intermedia vary between 161, 227 and 240 SWLI, and 178, 257 and 262 SWLI, respectively, in the Los Ríos, Chilo e and regional training sets). Figure 6 shows reconstructions of PMSE before and after the 1960 earthquake and highlights the similarity in model outputs at some sites ( Fig. 6A-D,G) and variance at others (Fig. 6E,F). Variance is in part due to a small number of diatom species having inconsistent elevation signatures across different modern sites and therefore, depending on the scale of the modern dataset used, differences in reconstructions can result (as described above at Llico). There are also problems associated with some species being unique to one particular sub-region, and therefore non-analogue situations arise in some models but not others, also leading to variation in reconstructions. For example, Luticola mutica and Tryblionella debilis constitute 13-30 and 7-13% of pre-and post-earthquake fossil assemblages, respectively, at Llico but these species are either absent or found in very low abundances in the sub-regional Chilo e modern dataset. Likewise, Fragilaria construens var. venter is abundant in fossil samples from Pucatrihue (up to 28%), and is equally abundant in the Chilo e training set (up to 35%), but constitutes <3% of any one sample in the Los Ríos training set. Figure 6 confirms that no model scale systematically produces higher or lower PMSEs, and also emphasizes that even where we see similarity in reconstructions between all scales, the models differ significantly when we assess the number of 'good' and 'close' modern analogues they provide. While local-scale models produce reconstructions with smaller error terms (45-64% of the difference between MHHW and MSL) than regional models (70-76%), they perform significantly worse when assessed using the MAT (Table 5). Regional models consistently provide the greatest number of 'good' modern analogues to fossil samples, with local scale models providing the fewest (Table 5; Fig. 6). We see an improvement in dissimilarity values with increasing spatial scale (Fig. 7). Figures 7A and B show the reduction in minDC values when using the regional or sub-regional training sets compared to limiting the training set to samples local to the core site only. Figure 7C shows the further advantage of using regional over sub-regional models.
Furthermore, variations in reconstructions at Chucal en between the two models both considered local (Table 5) highlight how local modern training sets should be used with caution. Chucal en is situated approximately equidistant between Quilo and Estero Guillingo, and both would be considered as local to the fossil core site (<3 km). However, if we relied upon local samples from only one of these sites, we could assume very different amounts of coseismic deformation occurred, due to differences in species found at both modern marshes as well as variations in elevation optima of principal diatom species between the sites. Eunotia paludosa and Pinnularia appendiculata together comprise 68-81% of pre-1960 diatom assemblages, and dominate the high marsh zone in the contemporary environment at Estero Guillingo, but are found in very low abundances (<5% total diatom valves counted) in only six modern samples from Quilo (Fig. 2). Consequently, elevation optima of these two species vary considerably between modern sites (e.g. by 85 SWLI units in the case of Pinnularia appendiculata between Quilo and Estero Guillingo), and cause significant variation in PMSE reconstructions. Furthermore, Fig. 4B highlights the sensitivity of species elevation optima to local site selection, which will cause inter-site variations in reconstructions. Elevation optima are far more consistent between both subregional and regional modern training sets (Fig. 4B). Of particular note at Chucal en is the influence of Navicula cincta on reconstructions, which constitutes up to 20% of post-earthquake diatom assemblages, and is seen to have a widely variable elevation optimum in different modern training sets. Figure 3 also confirms the inter-site differences between Estero Guillingo and Quilo, with the top-right cluster of Chilo e samples clearly forming two sub-groups, with different assemblages from similar elevations. Table 4. Ages of fossil samples determined by accelerator mass spectrometry radiocarbon dating. The radiocarbon chronologies of the Río Maullín and Chucal en sequences are reported in Cisternas et al. (2005) and Garrett et al. (2015), respectively.  (2017) Nevertheless, while this highlights that we should not rely solely on modern samples local to the fossil site, it is important to include them in regional models. When we analyse the similarity between fossil samples and their most similar modern samples, we see that for sites where we have samples from the local area in the regional modern training set, these provide most of the ten closest modern analogues (Fig. 8). In addition, at fossil sites where we have local modern samples, we see fossil samples have fewer 'good' or 'close' modern analogues when we exclude these local samples from the regional model. This is particularly apparent at Chaihuín, where the full regional training set provides all fossil samples with at least 'close' modern analogues; conversely, no fossil samples have even 'close' modern analogues if local samples are excluded from the regional training set (Table 5). This is a function of the unique species assemblages present at Chaihuín in modern samples from elevations above 240 SWLI (left-hand cluster in Fig. 3), which also dominate the fossil assemblages. In particular, key species (up to 30%) in fossil assemblages Rhopalodia gibberula and Eunotia praerupta are found in significantly higher abundances above 240 SWLI at Chaihuín than anywhere else (up to 34 and 10% of the modern assemblage, respectively, compared to <7 and 2%, respectively, elsewhere). This emphasizes that while regional models are best, they should include samples from the local area wherever possible.

Do Chilean tidal marshes underestimate coseismic deformation?
To assess the performance of models from different spatial scales, we compare our PMSE reconstructions with observations made 8 years after the earthquake of changes in pre-  Table 5. Land-level change associated with the 1960 earthquake, reconstructed using modern training sets from different spatial scales. In all cases we exclude modern samples <125 SWLI from transfer function models. Uplift is positive, subsidence is negative. We define 'good' and 'close' as the 5th and 20th percentiles, respectively, of dissimilarity values for the modern samples (Simpson, 2007;Watcham et al., 2013 (Plafker and Savage, 1970) not including possible local sediment compaction (Table 1). †Coseismic deformation is calculated as the difference between the elevation of pre-and post-earthquake samples (in metres above MSL) allowing for field elevation of each sample. The associated error term is the square root of the sum of the squared sample-specific standard errors of pre-and post-earthquake samples. Errors are AE2SD.  Figure 6. Diatom assemblages and palaeomarsh surface elevation reconstructions for the 1960 earthquake, showing comparison of results from local, sub-regional and regional scale transfer function models (WA-PLS component 2). Calibrated radiocarbon ages (in years AD) are shown where available (see Table 4 and references therein); Río Andali en and Río Tir ua were dated using caesium-137. Elevation optima of diatom species are classified based on weighted averaging of the regional modern training set. The 95% (2s) errors are shown for the regional model only to allow easier comparison of models on the left-hand graphs. minDC values are plotted relative to the threshold for a 'good' modern analogue (the 5th percentile and post-earthquake tide levels and vegetation growth limits (Plafker and Savage, 1970; Table 1). No model scale systematically agrees best with observations, and as numerical evaluation of models favours using regional-scale modern data, we focus here on comparing observations with reconstructions from the regional model (Fig. 9). Our transfer function produces reconstructions that fit with observations within errors at seven of the eight sites. It performs worst at two sites where there are no modern samples local to the fossil site included in the modern training set (Pucatrihue and Llico), for which there are three possible explanations. Firstly, the regional modern training set provides no fossil sample at either site with a 'good' modern analogue, and only provides 'close' modern analogues at Llico (Table 5; Fig. 6E,F). However, as the transfer function does well at similar sites with no local modern samples within the Bío Bío region (Fig. 6A,B), this suggests that diatom species distributions may be similar between the Bío Bío, Araucanía and Los Ríos regions and that some species may be unique to the Los Lagos region (e.g. Pinnularia obscura). Alternatively, the presence of the acidtolerant diatom P. obscura in fossil assemblages from Llico may indicate a unique pre-earthquake boggy environment, dominated by Sphagnum moss, which typically produces  (2017) organic acids (DeNicola, 2000), that is not captured in our modern sampling. Finally, it appears the dominance of two key diatom species account in part for the transfer function predicting less coseismic land-level change at Llico than observations suggest occurred. Denticula subtilis and Nitzschia fruticosa together comprise 14-22% of postearthquake fossil diatom assemblages at Llico, and our investigation of diatom elevation optima in the modern environment shows that WA-PLS component 2 significantly manipulates (increases) the optima of these species (Fig. 4A). As these species are only found post-earthquake, they would have the effect of increasing post-earthquake PMSE reconstructions, and therefore reducing the overall reconstructed coseismic land-level change.
The only other site where either D. subtilis or N. fruticosa are found in significant abundances in fossil assemblages is at  Río Andali en, where their presence could explain the discrepancy in the direction of land-level change between observations and predictions. D. subtilis is abundant both pre-and post-earthquake, but significantly more so preearthquake (up to 35% of the assemblage). This would make reconstructed pre-earthquake PMSEs abnormally high, resulting in apparent coseismic subsidence rather than the observed uplift.
The other site where our transfer function significantly underpredicts coseismic land-level change compared to observations is Río Lingue. However, the regional training set provides all fossil samples with at least a 'close' modern analogue, and provides 63% of samples with a 'good' modern analogue; instead the discrepancy could be due to inaccuracy in the original land-level change observation at this site, which Plafker and Savage (1970) note was made on unconsolidated sediments. Alternatively, the elevation of the fossil core site may be too low in the tidal frame given the problems the transfer function has predicting elevations below MSL (see above: 'Transfer function development' ;  Fig. 5). The transfer function reconstructs a pre-earthquake PMSE of $1.2 m above MSL, so if the observation of 1.6 m of subsidence is correct, then the location would be below MSL post-earthquake. We exclude samples with elevation optima below 125 SWLI (MSL ¼ 100 SWLI) from the modern training set, which therefore prevents the transfer function predicting elevations below MSL. However, even if these samples were included, we see a breakdown of species-elevation relationships in these samples, which would limit the ability of the transfer function to predict low elevations. Our reconstructions at Río Lingue could also be influenced by allochthonous diatom taxa, notably Paralia sulcata and Planothidium delicatulum var. delicatula, which account for $10-15% of total valves counted in pre-and post-earthquake samples, respectively. Due to Paralia sulcata's long-chained structure, its robust valves are known to be transported by tidal currents across a marsh and deposited far inland. This can therefore complicate our interpretations of the depositional environment, as its ecological preference is for habitats lower in the tidal frame than it is often found (Hemphill-Haley, 1995;Sawai et al., 2004;Dura et al., 2016). Likewise, the rapheless valve of Planothidium delicatulum is susceptible to transport up-marsh and may distort reconstructions in a similar way, resulting in higher reconstructed post-earthquake PMSEs than would be expected (Sawai, 2001;Sawai et al., 2016).
The influence of cosmopolitan diatom taxa could also help to explain why the transfer function underestimates coseismic land-level change compared to Plafker and Savage's (1970) observations (Table 5) at all sites except those in the Bío Bío region (Río Andali en, Río Tir ua). Several principal diatom species within the fossil assemblages with elevation optima in the modern environment below 150 SWLI have very broad tolerances of over 100 SWLI units (e.g. Planothidium delicatulum var. delicatula, Opephora pacifica, Opephora schulzii, Pseudostaurosira perminuta, Figs. 4A and 6), and therefore have the effect of disproportionately increasing PMSE reconstructions of low marsh environments (post-earthquake) compared to higher marsh environments (pre-earthquake). Nitzschia frustulum and Fragilaria construens var. venter also have similarly broad tolerances, and are significant drivers of the transfer function reconstructions at Llico and Pucatrihue, respectively.
Finally, there could be a delay in sediment deposition immediately post-earthquake, and therefore our estimates of coseismic subsidence may in fact include some post-seismic land uplift that has already occurred before sediment accumulation, resulting in underestimation of coseismic land-level change (Garrett et al., 2013;Dura et al., 2016). Garrett et al. (2013) found variability in sediment accumulation within the 2 years following the 2010 Maule earthquake, but observed little or no sediment accumulation at two out of four sites investigated. However, given that Plafker and Savage's (1970) observations were made in 1968, 8 years after the earthquake, these will also include post-seismic land-level change, and therefore any delay in post-seismic sedimentation will not be relevant when comparing with diatom-based observations. It therefore appears that the transfer function is underestimating coseismic deformation, and the dominance of allochthonous and cosmopolitan diatom species across Chilean tidal marshes provides a logical explanation.

Reconstructing PMSE in the absence of key data
To broaden the application of this investigation, we also evaluate model performance specifically at sites where we Figure 9. Comparison of observed (Plafker and Savage, 1970) and reconstructed landlevel change associated with the 1960 earthquake using the regional transfer function model (WA-PLS component 2). Copyright (2017) have no local modern diatom data or no observations of land-level change with which to compare our reconstructions. As Watcham et al. (2013) emphasized, these are the scenarios we most often face when attempting relative sea-or land-level reconstruction. At the three sites where we have a local modern training set, we also re-ran the regional transfer function excluding local samples (Río Lingue, Chaihuín, Chucal en; Table 5). The difference between reconstructions using the regional model and regional model excluding local samples is a maximum of 0.26 m at Chaihuín, and <0.01 m at Río Lingue and Chucal en. While excluding local samples from the model provided fossil samples with fewer 'good' or 'close' modern analogues than the full regional model, and we also concluded from Fig. 8 that modern training sets should include samples from the local area wherever possible, the variation in reconstructions is well within the error terms of the reconstructions. Therefore, we advocate that while it is still possible to reconstruct PMSE in the absence of local modern data, we must exercise caution in doing so. Where the collection of modern samples is simply not possible due to resource constraints or a restricted range of contemporary environments, we should investigate the elevation optima of principal diatom species, and only have confidence in the reconstructions where the transfer function does not extensively manipulate the data. In such situations, the key question is whether sub-regional or regional modern training sets are more appropriate. We therefore looked at sites with no local modern samples (Río Andali en, Río Tir ua, Pucatrihue, Llico and Río Maullín), and found the closest modern analogues to fossil samples do not necessarily come from the closest site, or even closest sub-region (Fig. 8). If using a sub-regional model, logically the sub-region closest to the fossil site would be chosen, but if we used the subregional model closest to the fossil sites at Llico (Chilo e) and Pucatrihue (Los Ríos), we would lose over 70% of the closest modern analogues. This therefore supports our earlier interpretations that regional models are superior to sub-regional or local ones.
The other key question besides what to do in the absence of local modern data is how far back in time does any part of the local marsh no longer provide a modern analogue? To address this we use the Chucal en fossil record that covers the last four great earthquakes that have been inferred for the 1960 region (Table 6; Fig. 10). Throughout the entire 1000year sequence, the regional model performs best in terms of providing fossil samples with the best modern analogues. Even for the oldest event, dated to AD 1070-1220, the regional modern training set provides 69% of fossil samples with at least a 'close' modern analogue, and 31% with a 'good' modern analogue. By contrast, the sub-regional Chilo e model does not provide fossil samples with any 'good' modern analogues for any of the four palaeoearthquakes. We do not observe increased variation between models or a decrease in 'good'/'close' modern analogues as we go further back in time. For three events, we also see there are significant differences between the two local models (Quilo and Estero Guillingo), because of variation in both species composition and elevation optima of principal species between the two local modern training sets (as discussed above: 'Assessing the optimum spatial scale of the modern training set'). This further highlights the potential danger of basing reconstructions of coseismic deformation on local samples. It is therefore important not only to have as many modern samples from the local area as possible, but also not to rely solely on modern samples from the local or sub-regional area when reconstructing deformation beyond recent events for which there are observations. Likewise, it is essential to understand the diatom ecology and explore the potential impact of key species on reconstructions.

Conclusions
This paper offers advances in improving understanding of (i) modern diatom distributions in south-central Chile, in the  (2017) region of the 1960 earthquake and (ii) the diatom-based transfer function methods that we will use in the future to provide quantitative estimates of deformation that occurred in pre-1960 earthquakes preserved in the geological record. Ultimately, it is not possible to constrain earthquake sizes and spatial patterns of land movement from tidal marsh records without this improved knowledge.

Land-level change in 1960
We have demonstrated the viability of our transfer function methods by using them to reconstruct land-level change associated with the giant 1960 earthquake at five new sites, confirmed reconstructions at an additional three sites and provided reconstructions that match observations. We have extensively explored transfer function model outputs, comparing reconstructions derived from different scales of modern data. We have assessed transfer function performance through numerical evaluation of transfer function models and using the modern analogue technique, but as this can fail to discriminate between acceptable and unacceptable reconstructions even when there are good modern analogues (e.g. Woodroffe, 2009), we also advocate investigating elevation signatures of principal diatom species to assess the reliability of reconstructions. Figure 11 provides a summary of land-level change in 1960 reconstructed using our regional diatom transfer function. The general profile of vertical deformation from north to south agrees well with the observations of Plafker and Savage (1970), although reconstructions of coseismic subsidence are up to $90 cm less than observations. This is partly due to the dominance of allochthonous diatom species and the inability of the transfer function to predict elevations close to, and below, MSL. We also demonstrate that obtaining accurate reconstructions can be problematic when species whose elevation optima are sensitive to WA-PLS regression make significant contributions to the fossil dataset. Transfer function reconstructions from Chilean tidal marshes must therefore be viewed as minimum estimates of coseismic land-level change, and it is essential to investigate the elevation optima of principal diatom species in the modern environment and how these impact upon reconstructions.
We have also made advances in reducing the error terms on the reconstructions compared to previous studies by using Figure 10. Relative sea-level change over 1000 years at Chucal en (CH11/1), showing differences in reconstructions and 'good' modern analogues between regional, sub-regional (Chilo e), and local (Quilo) transfer function models. Elevation optima of diatom species are classified based on weighted averaging of the modern training set. minDC values are plotted relative to the threshold for a 'good' modern analogue (the 5th percentile). The 95% (2s) error bars are only shown for reconstructions using the regional model to allow easier comparison between models. Earthquake horizons A, B, C and D are dated to AD 1960AD , 1575AD , 1270AD -1450AD and 1070AD -1220, respectively, as detailed by Garrett et al. (2015). The minerogenic unit overlying soil C is not interpreted as being a tsunami sand sheet (Garrett et al., 2015), and therefore reconstructions are shown throughout this unit (but not shown through inferred tsunami deposits A, B and D). an expanded modern training set. For those sites where reconstructions have been previously published at Río Andali en and Río Tir ua, we have reduced transfer function errors by 42 and 40%, respectively (cf. Garrett et al., 2013). At Chucal en, where the record extends beyond 1960 to include three earlier palaeoearthquakes, our new modern diatom data help confirm previously published estimates of coseismic deformation for all four events, although fossil samples have better modern analogues and consequently transfer function errors are reduced by between 28 and 40% (cf. Garrett et al., 2015). While our reconstructions are sufficiently precise to provide important constraints on coseismic land-level change, error terms could still be considered high compared to tidal ranges and compared to the magnitude of reconstructed coseismic land-level change. Consequently, we make the following recommendations for future research.

Lessons from Chile and recommendations
We present an expanded modern diatom data set for southcentral Chile showing vertical zonation of species, particularly above MHHW. We model the contemporary relationship between species and elevation to allow us to quantify past changes in marsh surface elevation, and investigate the effects of using modern training sets from local, sub-regional and regional scales on reconstructions and model performance. We conclude: Regional models are superior to local models, both in terms of agreement with observations and in providing fossil samples with the closest modern analogues.
Modern training sets should include samples from the local area where possible, as local samples provided the majority of the ten closest modern analogues.
In the absence of local modern samples, regional modern training sets should be used in preference to sub-regional sets, as the closest modern analogues do not always come from the closest sub-region. Modern training sets should include samples from as broad an area as possible.
Our findings that regional models out-perform subregional and local scale models agrees with previous studies from North America (Nelson et al., 2008;Watcham et al., 2013;Sawai et al., 2016) and the UK (Gehrels et al., 2001;Horton and Edwards, 2005), which have also highlighted the need for large modern training sets that include samples from the full environmental gradient. However, it does appear that our Chilean transfer function does not perform as well as its regional counterpart in Alaska , in terms of providing fossil samples with as many 'good' modern analogues and producing error terms that are a smaller percentage of the tidal range. The reason for the former may simply be a function of modern training set size, as the regional Alaskan training set consists of 255 samples, compared to 176 in this paper. However, the latter appears to be more of a genuine problem, with diatoms not being as good at reconstructing elevations below MHHW in Chile due to many dominant species having broad elevation tolerances between 0 and 200 SWLI units (MHHW).
In future, when reconstructing deformation associated with palaeoearthquakes in Chile we should therefore (i) collect new modern samples to further expand the regional modern training set, and (ii) target fossil sites higher in the tidal frame when we expect coseismic subsidence has occurred. The latter should increase the likelihood of reconstructing elevations where diatom distributions are more tightly defined than they are around present MSL, and therefore reduce the influence of species with broad environmental tolerances that reduce the precision of transfer function reconstructions. In areas where low-elevation diatom taxa are under-represented on modern marshes, or where taphonomic processes occur at lower elevations, problems of transfer function under-prediction may remain. The final recommendation is therefore to ensure careful examination of diatom elevation optima and tolerances to highlight such instances where reconstructions can only provide minimum estimates of land-level change. All these recommendations are also pertinent when reconstructing seismic land-level change on subduction zone coasts worldwide.
Acknowledgements. This research was funded by the Natural Environment Research Council (New Investigator Award NE/ K000446/1 to E.H. and E.G.) and FONDECYT (Project No. 1150321 to M.C.). Radiocarbon dating support was provided by the Natural Environment Research Council Radiocarbon Facility (grants 1707.0413 and 1795.0414), and we gratefully acknowledge Steve Moreton for providing the radiocarbon dates. We thank Bill Austin, Sarah Woodroffe and Caroline Taylor-Garrett for help in the field, and thank Ian Shennan for providing helpful comments to improve an earlier draft of this paper. We appreciate the constructive reviews of Robin Edwards and an anonymous reviewer, which greatly improved the original manuscript. This paper is a contribution to IGCP Project 639. Table S1. Stratigraphy of the four newly cored fossil sites. Cores in bold are analysed in this paper. Sediments were described and then subsequently interpreted as being of seismic origin following the criteria of Nelson et al. (1996). We employed the criteria of Peters and Jaffe (2010) to confirm sand layers as being indicative of tsunami deposition. For the stratigraphy of previously published sites, see Garrett et al. (2013) and Cisternas et al. (2005). Table S2. Key diatom species and their salinity and ecological preferences. All taxa listed exceed 5% of total  (2017) diatom valves counted in at least one sample, with those exceeding 20% of at least one modern sample in bold. Some species do not exceed 5% in both modern and fossil samples, as indicated by the presence in brackets. Ecological information and salinity preferences from Hassan et al. (2009), Hemphill-Haley (1993, Van Dam et al. (1994) and Vos andDe Wolf (1988, 1993). Figure S1. Modern tidal marsh characteristics, including sampling transects, marsh zonation based on vegetation, dominant plant species and dominant diatom taxa exceeding 15% of total diatom valves counted in at least one sample. Species are ordered by elevation optima, with dark and light blue reflecting elevation optima below and above MHHW, respectively.