Lithology and microfossil biostratigraphy beneath the marshes of a central Oregon estuary limit geophysical models of Cascadia megathrust rupture during successive earthquakes by ruling out >0.5 m of coseismic coastal subsidence for the past 2000 yr. Although the stratigraphy in cores and outcrops includes as many as 12 peat-mud contacts, like those commonly inferred to record subsidence during megathrust earthquakes, mapping, qualitative diatom analysis, foraminiferal transfer function analysis, and 14C dating of the contacts failed to confirm that any contacts formed through subsidence during great earthquakes. Based on the youngest peat-mud contact’s distinctness, >400 m distribution, ∼0.6 m depth, and overlying probable tsunami deposit, we attribute it to the great 1700 CE Cascadia earthquake and(or) its accompanying tsunami. Minimal changes in diatom assemblages from below the contact to above its probable tsunami deposit suggest that the lower of several foraminiferal transfer function reconstructions of coseismic subsidence across the contact (0.1–0.5 m) is most accurate. The more limited stratigraphic extent and minimal changes in lithology, foraminifera, and(or) diatom assemblages across the other 11 peat-mud contacts are insufficient to distinguish them from contacts formed through small, gradual, or localized changes in tide levels during river floods, storm surges, and gradual sea-level rise. Although no data preclude any contacts from being synchronous with a megathrust earthquake, the evidence is equally consistent with all contacts recording relative sea-level changes below the ∼0.5 m detection threshold for distinguishing coseismic from nonseismic changes.

The series of great (moment magnitude >8) earthquakes early in the twenty-first century has resulted in greater appreciation for the variability of megathrust earthquake ruptures at subduction zones (Wang, 2007; Melnick et al., 2012; Wang and Tréhu, 2016; Bilek and Lay, 2018; Wang et al., 2018). Because such variability complicates the local as well as ocean-wide earthquake and tsunami hazard forecasts used to direct hazard mitigation, reconstructing the history of the greatest ruptures and their accompanying destructive tsunamis remains fundamental to hazard assessment (Mueller et al., 2015; Wirth and Frankel, 2019). Once used primarily to estimate the average recurrence of great earthquakes for entire subduction zones, the chief benefit of recently developed earthquake and tsunami histories is to limit increasingly complex models of megathrust rupture to what has happened in the past (e.g., Witter et al., 2012; Nelson, 2013; Wang et al., 2013; Moernaut et al., 2014; Shennan et al., 2016; Gao et al., 2018; Wirth and Frankel, 2019). The most valuable histories—particularly in subduction zones that lack long historical records—include reconstructions that extend models based on instrumental measurements back in time through multiple cycles of great earthquakes (e.g., Ely et al., 2014; Shennan et al., 2014; Garrett et al., 2016; Hayward et al., 2015; Meltzner et al., 2015; Wesson et al., 2015; Milker et al., 2016; Pinegina et al., 2020). Coast-based histories are too far landward to rule out competing models of megathrust rupture (Wang et al., 2013; Wang and Tréhu, 2016). However, such histories limit models by showing differences in the amount of overriding plate deformation during successive earthquakes (or the inundation extent of their tsunamis), either (1) along different segments of a subduction zone at about the same time (e.g., Leonard et al., 2004; Van Daele et al., 2015; Wang et al., 2013; Shennan et al., 2014; Kemp et al., 2018; Padgett, 2019) or (2) over time at a site (Nelson, 2013; Dura et al., 2016a; Sawai et al., 2004; Cisternas et al., 2005; Enkin et al., 2013; Briggs et al., 2014; Clark et al., 2015; Shennan et al., 2016; Dura et al., 2017; Moernaut et al., 2018; Hong, 2019). Here, we describe evidence at a coastal site in the central Cascadia subduction zone where lithostratigraphy and biostratigraphy set limits on models of successive earthquake ruptures by ruling out substantial coastal subsidence for the past 2000 yr.

Following studies of coastal deformation during magnitude 9 subduction-zone earthquakes in Alaska (Plafker, 1969; Ovenshine et al., 1976; Bartsch-Winkler et al., 1983) and Chile (Wright and Mella, 1963; Plafker and Savage, 1970; Plafker, 1972), the interpretation of coastal wetland stratigraphy along the coasts of Washington, Oregon, and northern California as an archive of regional vertical deformation during great earthquakes has helped to end debate about whether or not Cascadia’s subduction-zone megathrust slips smoothly or has been locked for hundreds of years and is storing strain to be released in a future great earthquake (Savage et al., 1981; Heaton and Kanamori, 1984; Adams, 1984; Heaton and Hartzell, 1987; Atwater, 1987; West and McCrumb, 1988; Darienzo et al., 1994; Nelson and Personius, 1996). A key tenant of early studies was that tidal wetland stratigraphy of Cascadia—where the coast experiences successive cycles of megathrust overriding-plate deformation—differed from the stratigraphy beneath similar temperate wetlands along passive-margin coasts. Based on initial studies (e.g., Atwater, 1987, 1992; Darienzo and Peterson, 1990; Nelson, 1992a; Clarke and Carver, 1992), the sharp stratigraphic contacts between thin beds of peaty sediment of former marshes and swamps overlain by much thicker beds of muddy (rarely sandy) tidal-flat sediment were inferred to record the jerky rise of late Holocene relative sea level (RSL) punctuated by sudden subsidence during successive great megathrust earthquakes (Atwater et al., 1995; Nelson and Personius, 1996). This tidal stratigraphy of interbedded lithologies reflecting jerky RSL rise was contrasted with the 1- to 4-m-thick sections of largely peaty wetland sediment common on temperate North American coasts, which were interpreted as the product of gradual late Holocene sea-level rise (e.g., Bloom and Stuiver, 1963; Redfield, 1972). Although some cautioned that vertical tectonic deformation is only one of many factors that influence tidal sedimentation at Cascadia (Darienzo and Peterson, 1990; Nelson, 1992b; Long and Shennan, 1994; Nelson et al., 1996b; Allen, 2000), the model of jerky late Holocene RSL rise remained the basis for interpretations of repeated subsidence of tidal wetlands during as many as 12 great earthquakes at tens of sites along the subduction zone (Atwater et al., 1995; Nelson and Personius, 1996; Clague, 1997; Shennan et al., 1998; Kelsey et al., 2002; Witter et al., 2003; Nelson et al., 2004; Schlichting and Peterson, 2006; McCalpin and Carver, 2009; Valentine et al., 2012; Graehl et al., 2014; Hutchinson and Clague, 2017; Hong, 2019; Padgett, 2019; Nelson et al., 2020). Similar assumptions were used to infer that a record of megathrust earthquake deformation is preserved in tidal sequences on other subduction-zone coasts (Nelson, 2013; Dura et al., 2016a; Shennan et al., 2016).

However, the largely successful application of the jerky RSL rise model at Cascadia had two unintended consequences. The first was that it helped to obscure significant along-strike differences in tidal stratigraphy—likely reflecting differences in RSL and(or) earthquake history—along the subduction zone (Nelson and Personius, 1996; Nelson, 1992b). The second was that with uncertainties in radiocarbon dating of many decades to centuries (e.g., Nelson, 1992a; Graehl et al., 2014; Hutchinson and Clague, 2017), it fostered correlation of sharp stratigraphic contacts for many hundreds of kilometers along the subduction zone. The latter, in turn, led investigators to infer—or at least prevented them from discounting (e.g., Atwater et al., 1991; Nelson et al., 1995)—an earthquake history of primarily long ruptures during earthquakes near magnitude 9. It remains unresolved for most coastal sites whether such a history of mostly giant earthquakes is the result of an actual difference in rupture history, unlike that of other subduction zones (e.g., Wang et al., 2013; Wang and Tréhu, 2016; Bilek and Lay, 2018), or a lack of preservation of coastal evidence for ruptures of a few hundred kilometers or less (Nelson et al., 2006; Shennan et al., 2016; Hutchinson and Clague, 2017).

In this paper, we describe the stratigraphy in cores and outcrops at the Siuslaw River estuary of central Oregon that is much like those commonly inferred to record a series of megathrust earthquakes on the temperate coasts of this and other subduction zones. The sequence includes 9–12 peat-mud contacts that potentially record more earthquakes in the past 2000 yr than at any of the tens of tidal sites to the north and south (Figs. 1 and 2). If attributed to earthquakes, such a stratigraphy might steer debate about the frequency and coastal extent of past great earthquakes at Cascadia (Nelson et al., 2006; Frankel, 2011; Goldfinger et al., 2012, 2016; Atwater et al., 2014; Hutchinson and Clague, 2017). Instead, our mapping of stratigraphic contacts beneath tidal marshes near the river, lithologic descriptions of cores and outcrops, qualitative diatom analysis, quantitative foraminiferal analysis using a Bayesian transfer function, and 14C dating of most of the contacts failed to confirm that any of the contacts formed through sudden subsidence during great earthquakes. The failure, however, constrains models of megathrust rupture to ∼0.5 m or less of coseismic subsidence along this part of the Oregon coast for the past 2000 yr. The failure also shows the utility of the criteria of Nelson et al. (1996a) and Shennan et al. (2016) in identifying earthquake contacts, and it further illustrates how thresholds (e.g., McCalpin and Nelson, 2009) for the creation and preservation of earthquake contacts limit their identification at the Siuslaw River estuary and, by analogy, at similar sites elsewhere.

Detailed location maps, figures showing additional stratigraphy, tables of detailed data, standard methods (such as measurement of elevations and tide levels), and summaries of the tidal marsh setting of our study area and the two previous investigations of Siuslaw River stratigraphy (e.g., Fig. 3; Nelson, 1992b; Briggs, 1994) appear in the Supplemental Material1 for this paper (Parts 1–3).

Mapping Potential Earthquake Contacts

To reevaluate the tidal stratigraphy of the Siuslaw River estuary—and the alternative interpretations of it (Figs. 2 and 3; Supplemental Material Part 3)—in the context of recent studies of great earthquake and tsunami stratigraphy at Cascadia, in 2007–2009 we examined the interbedded sequences of peaty and muddy sediment of northern and eastern Cox Island in greater detail (Figs. 47). In particular, with more detailed lithologic descriptions of many more cores, improved methods of microfossil analysis, and more precise 14C dating, we sought to identify sudden changes in RSL that might correlate with evidence of earthquakes and tsunamis at more recently studied coastal sites and earthquake-generated turbidites offshore (e.g., Graehl et al., 2014; Goldfinger et al., 2016; Milker et al., 2016; Nelson et al., 2020). We examined 61 gouge cores (100-cm-long segments, 25 mm diameter), describing 38 of them in the field along with stratigraphic sections at outcrops on the north side of Cox Island and 200 m upriver from its northeast corner using the Troels-Smith (1955) system for describing organic-rich sediment (methods described by Nelson, 2015; see Figs. 2, 4, 5, and 6; Fig. S5 [see footnote 1]). Four vibracores (2- to 5-m-long continuous cores, 70 mm diameter) were collected for 14C dating and microfossil analysis at two sites (S and K, Figs. 4 and 5).

The vibracores with the most distinct stratigraphy at each of the two sites were split, photographed, and wrapped for transport the day after collection. In the laboratory, we described their lithostratigraphy in detail and noted other sediment characteristics and sediment color with Munsell color charts (Fig. 5; e.g., Nelson, 2015). At site S, the first vibracore (Sa) was sampled for diatoms and foraminifera; a second vibracore (Sb), 4 m from the first, was used for radiocarbon and later foraminiferal sampling. Except when discussing the later foraminiferal samples from core Sb (Table S2), we combined the stratigraphy and samples from the two vibracores at each site, referring to them as cores S and K. Depths of lithologies and contacts in cores S and K approximately corrected for 20%–25% compaction are shown in Figures 5 and 6.

On Cox Island, 38% ± 18% (error = 1σ) of units in 13 representative cores along three core transects and in two stratigraphic sections at outcrops consisted of peaty beds (Table 1; Fig. 4). Freshly exposed peat units (estimated organic component by volume >50%) typically showed 5YR–7.5YR color hues, muddy peat commonly had 10YR hues, and the hues of muddy units ranged from 2.5Y to 5Y, commonly with lighter color values than darker peaty units. By comparing lithologic descriptions and photomosaics of vibracores and some gouge cores, we traced the nine most continuous of 12–15 peaty beds (0–4 m depth) for 100–400 m (Fig. 7). Correlation of intervening less peaty beds over >50–100 m was less certain. Twice as many (42% vs. 23%) upper contacts on the nine primary peaty beds were sharp (<3 mm) compared with lower bed contacts. However, only half of lower contacts were gradational enough to suggest gradual marsh emergence, and a third of upper contacts were more consistent with gradual rather than sudden submergence. We labeled the nine primary peat-mud contacts at the tops of peaty beds A through I (Figs. 5, 6, and 7; see discussion later herein). Three distinct but less extensive peat-mud contacts, which were used with 14C ages to correlate cores S and K with outcrop 1, led to subdividing contacts D, E, and F into contacts Da, Db, Ea, Eb, Fa, and Fb (Figs. 5 and 6). Contacts B, Ea, Fa, and H in the vibracores were not identified in the outcrop.

Modeling Contact Ages

We estimated the times contacts formed using accelerator mass spectrometer (AMS) 14C ages on samples of plant fragments from above and below contacts in cores S and K, and in blocks of sediment spanning the contacts cut from outcrop 1 (Figs. 4, 6, and 8; Table 2). Most of the 60 samples were selected by washing 3- to 5-mm-thick vertical slices of sediment on a 1 mm sieve under a binocular microscope (6–50×; methods of Kemp et al., 2013).

We used OxCal stratigraphic ordering software (Bronk Ramsey, 2008, 2009) to develop a series of age models for the 12 contacts (nine primary contacts of Figure 7 and three less extensive contacts identified in cores S and K and at outcrop 1 as shown on Figure 6; Table 2). Initial modeling consisted of outlier analyses (methods of Bronk Ramsey, 2009) starting with all ages, most grouped into OxCal phases (groups consisting of unordered samples) above and below contacts (OxCal code for selected models in Supplementary Material Part 5 [footnote 1]). In the series of age models, we then successively eliminated ages that were obvious outliers or that we interpreted to be less accurate minimum or maximum estimates of the times contacts formed (e.g., Milker et al., 2016; Witter et al., 2019; Nelson et al., 2020). Our inferred closest maximum and minimum ages for each contact are marked in bold on Table 2. For our final age model, we used an OxCal sequence (nonoutlier) analysis model with only the closest (youngest) maximum age and(or) closest (oldest) minimum age for each contact (ages marked by asterisks on Table 2). As it is unlikely that the different types of dated materials were from the same age population, the closest ages better restricted modeled ages for each contact than did averages of similar ages (e.g., Johnstone et al., 2019; Streig et al., 2020).

We based our interpretations of the closest maximum and minimum ages (discussed for each contact below) on the type of plant macrofossil, its orientation, degree of decay and abrasion, host stratigraphic unit lithology, its stratigraphic context relative to adjacent plant macrofossils and to upper and lower units, and—most importantly—its calibrated 14C age relative to the ages of samples above and below the sample of interest. As elsewhere in Cascadia coastal sequences (e.g., Nelson et al., 2006; Hutchinson and Clague, 2017), most of our ages were on detrital materials, which are older than the times at which adjacent contacts formed. The relative ages of rhizomes (belowground stems) of low and middle marsh herbs are more difficult to interpret than ages on aboveground plant parts. Usually growth-position rhizomes, especially those of Triglochin maritima with the bases of its decay-resistant leaves still attached, provide unambiguous minimum ages for underlying contacts. Rarely, we inferred from the sequence of ages on adjacent samples that the rhizomes of plants younger than contacts grew down into the peaty unit just below a contact (samples OS-138531, OS-62145, OS-62219, OS-66499, OS-144809; Table 2) and, therefore, provide minimum ages for an overlying contact.

We also dated materials in 11 samples from six gouge cores (7, 8, 9, 11, 12, 19) collected in 1987 (Figs. 2 and 4; Figs. S2, S3, and S1H; Table 2). Although one far-too-young, outlier age (655 ± 15 14C yr B.P.) indicated needles dragged down from a higher level in the core, the other ages were consistent with the recent ages of Table 2 (5 of the 11 ages are shown in fig. 2 of Nelson, 1992b).

Microfossil-Based Assessments of Environmental Change across Contacts

Over the past two decades, the use of changes in fossil foraminiferal and diatom assemblages to stratigraphically identify great earthquakes at Cascadia has shifted from using mostly qualitative and limited quantitative comparisons of assemblages to estimate amounts and rates of RSL rise across peat-mud contacts (e.g., Nelson et al., 1996b, 1998; Atwater and Hemphill-Haley, 1997; Shennan et al., 1998; Kelsey et al., 2002; Witter et al., 2003; Hawkes et al., 2005; Graehl et al., 2014; Hemphill-Haley et al., 2019) to transfer function methods that produce sample-specific errors (Dura et al., 2016b; Guilbault et al., 1995, 1996; Nelson et al., 2008; Hawkes et al., 2011; Milker et al., 2016; Shennan et al., 2016; Horton et al., 2017). Transfer functions use the relations among modern assemblages and their respective elevations in modern tidal environments as analogs to hindcast past tidal elevations from fossil assemblages in stratigraphic sequences (Kemp and Telford, 2015). The most recent developments are Bayesian transfer functions (Cahill et al., 2016; Kemp et al., 2018; Hong, 2019; Padgett, 2019; Nelson et al., 2020), which allow species response curves to deviate from a predefined form (commonly unimodal) and may incorporate prior information about sampled sediment (i.e., stratigraphy, lithology, paleoecologic information from other types of fossils) to help constrain reconstructions of past RSL change.

Foraminiferal Analysis

At the Siuslaw River estuary, Kemp et al. (2018) used the original foraminiferal data of Hawkes et al. (2011) (see Table S2, 15 samples near contact A) with a new Bayesian transfer function to estimate the amount of rapid submergence (inferred to be the result of coseismic subsidence) across contact A in core S on Cox Island (Figs. 4, 5, and 9). Kemp et al.’s (2018) much larger data set (393 samples from 19 sites) than used to develop previous, non-Bayesian transfer functions (e.g., Hawkes et al., 2010; Engelhart et al., 2013a, 2013b; Milker et al., 2015b, 2016) included modern assemblages that are better analogs for fossil assemblages than those of earlier studies. With the new foraminiferal data reported here (22 samples), we used the same transfer function (informed West Coast function of Kemp et al., 2018) to reconstruct the amount of subsidence across contacts A, Db, Fa, and I in core S, and contacts C and Fb at outcrop 1 (Figs. 5, 6, 9, and 10). The 48 samples across other contacts (B, Da, Ea, Eb, G, and H) were barren or contained too few foraminifera for reliable subsidence estimates (Table S2).

All 85 samples of foraminifera were refrigerated, prepared, and counted using standard methods (e.g., Scott and Hermelin, 1993; de Rijk, 1995; Kemp et al., 2009; Engelhart et al., 2013b; Milker et al., 2015a). Core S was sampled across contacts A, Db, Fa, and I eight months after collection (2008). However, contacts Ea, Eb, G, and H were not sampled until a decade later; the later samples came from the second (refrigerated) vibracore Sb. Similarly, samples across contacts C, Eb, Fb, and G were from (refrigerated) blocks cut from outcrop 1 but were not sampled until 2019. Although 56% of the foraminifera samples were from sediment refrigerated for a decade prior to sampling, and 40% of the 85 samples were barren of foraminifera, we observed no tendency for the samples analyzed in 2019 to contain lower concentrations of foraminifera than those analyzed in 2008. Although some samples from core Sb were barren, other samples from similar lithologies in the same sections of core Sb had concentrations as high as those in adjacent core Sa, sampled in 2008 (Table S2). Ten genera and species of foraminifera were identified using the taxonomic illustrations and descriptions in Horton and Edwards (2006), Hawkes et al. (2010), Wright et al. (2011), and Milker et al. (2015a) (see Table S2).

To make our transfer function reconstructions of RSL change consistent with those of Kemp et al. (2018), we followed their procedures. We standardized our taxonomy, which differs slightly from the taxonomy for contact A of Hawkes et al. (2011) (see Table S2), by renaming Trochamminita irregularis to Trochamminita sp., and by combining all species of Haplophragmoides (Haplophragmoides maniliensis and Haplophragmoides wilberti) and calcareous species into single groups, respectively. Similarly, we excluded assemblages with <30 foraminifera from the reconstructions (Table S2; e.g., Hawkes et al., 2011; Kemp et al., 2018) because they may not be in situ assemblages, or they may have undergone significant taphonomic change and thus are likely to be unrepresentative of the environment at the time they were deposited. To check that our sample assemblages had good modern analogs in the Kemp et al. (2018) data set, we used the same modern analog evaluation technique: All but four of our fossil samples containing >30 foraminifera (at 281, 283, 324, 326 cm depths in core S; Table S2) met a 10% dissimilarity threshold in pairwise comparisons. The deeper two of those samples contained >96% Acostata mariae, a tidal-flat species not well represented in the Kemp et al. (2018) data set. Using “SWLI,” a standard water level index that allows comparison among sites with differing tidal ranges (e.g., Horton and Edwards, 2006; Kemp and Telford, 2015), we equated mean higher high water (MHHW) with 200 SWLI and mean tide level (MTL) with 100 SWLI, which near Cox Island are 1.19 m (2.30 m in North American Vertical Datum of 1988 [NAVD88]) and 0.02 m MTL (1.12 m NAVD88), respectively.

A key aspect of Kemp et al.’s (2018) new Bayesian transfer function is that it includes prior information about sample lithology (Figs. 5, 6, 9, and 10; Table S2; Cahill et al., 2016). Following general inferences about the elevational range of tidal sediment used in almost all studies of tidal stratigraphy at Cascadia, clastic-dominated samples typical of tidal flats or low marshes are assumed to have accumulated between local mean low water (18.1 SWLI or −0.98 m MTL at Siuslaw River) and MHHW (200 SWLI or 1.19 m MTL). Alternatively, organic-rich sediment, which commonly reflects middle and high tidal marsh settings, is assumed to have accreted above local mean high water (182 SWLI or 0.98 m MTL). The upper bound of the latter is the highest occurrence of foraminifera in the Kemp et al. (2018) data set (252 SWLI). Inclusion of the two lithologic priors in the transfer function analysis influenced reconstructed sample elevations by specifying that they were more likely to fall within the range of the assigned lithologic prior. These priors overlap and are conservative, in that they allowed the function to reconstruct RSL changes reflecting either submergence or emergence (Kemp et al., 2018). As discussed below for four samples of muddy peat above contact A, because the lithology of the samples suggests that they were deposited within the elevational range of the overlap between mean high water (MHW) and MHHW, assignment of the lower lithologic prior versus the higher prior to the samples gives differing results for subsidence across contact A (Figs. 9A and 9B). Because the lithology of the other samples does not suggest that they were deposited within the range of the overlap of the two priors, the uncertainty in which of the two lithologic prior groups to assign to samples does not apply to our other subsidence reconstructions (Figs. 9 and 10; Table S2).

Our reconstructions of subsidence across peat-mud contacts calculated with the Bayesian transfer function were so low that three of five 1σ errors and all 2σ errors on the reconstructions included negative values (Figs. 9 and 10; Table S2), which imply that the contacts could mark either uplift or subsidence. However, as widely assumed for decades (Nelson et al., 1996a) and recently shown by Horton et al. (2018) for Great Britain, peat-to-mud contacts in Holocene coastal sequences are far more likely to record submergence rather than emergence. Although the Kemp et al. (2018) Bayesian transfer function does not incorporate this assumption into subsidence calculations (Table S2), we inferred that only the positive intervals of our reconstruction errors were accurate (indicating submergence). This reduced the range of our reconstruction errors to less than those of most other similar studies (e.g., Kemp et al., 2018; Padgett, 2019).

Diatom Analysis

Diatoms in core S showed a more complete history of paleoecologic change than foraminifera because diatom samples came from longer sections of the core (Fig. 11). In 2008, we collected 136 4- to 7-mm-thick samples at 1 to 4 cm intervals above and below contacts, except contact Fb (Fig. 11; Tables S3 and S4). Diatom slides were prepared and counted using standard methods (e.g., Sawai et al., 2002; Sawai and Nagumo, 2003). About 250 diatom valves were identified in each sample under an oil-immersion microscope at 600× magnification (Table S4), including 258 taxa in 70 genera (Table S3). Fragments containing more than half a valve were included in the counts. Diatom abundance is shown as a percentage of the total number of diatom valves counted, with only taxa that exceeded 5% of valves in more than five samples used for paleoecological interpretation (Table S4). Because these criteria yielded 53 taxa, and meaningful summaries of changes in tidal diatom assemblages are complex (Dura et al., 2016b; Shennan et al., 2016), on Figure 11 we show only the 28 taxa that exceeded 10% of valves in six or more samples.

As have fossil diatom studies with similar objectives at similar sites (e.g., Sawai et al., 2004; Nelson et al., 2008; Shennan et al., 2016; Watcham et al., 2013; Dura et al., 2016b), we attempted to use a diatom transfer function to reconstruct past RSL for fossil diatom samples in core S, as we did with the foraminifera. In developing our diatom transfer function, we used the modern diatom data set of Sawai et al. (2016), which included 175 diatom assemblages from nine tidal marshes in Oregon and a tidal marsh in southwest Washington. Following Kemp et al.’s (2018) foraminiferal transfer function analysis, we used the SWLI index to standardize our elevations relative to local MHHW (SWLI = 200) and mean lower low water (MLLW) (SWLI = 18.1) to account for the wide variation in tidal range for sampled marshes in the modern data set.

Our diatom transfer function analysis followed routine procedures widely used in such analyses (Kemp and Telford, 2015). For our analysis, we applied a weighted averaging–partial least squares transfer function (Fig. S6). To improve the performance of the transfer function, we removed one of the assemblages in the modern data set. In the fossil assemblages, we excluded 81 of the 258 taxa not present in the modern data set of Sawai et al. (2016), as well as taxa for which maximum percentages were <2%. We also removed planktonic taxa (Aulacoseira, Skeletonema, Thalassionema, and Thalassiosira). The distributions of most planktonic species are not controlled by elevation but by environmental variables, such as salinity and pH. Their frustules and valves are easily transported by riverine and tidal currents, which are not dependent on sampling elevation. Although taxa of Melosira are sometimes classified as planktonic, we retained this group in our analysis, following Sawai et al. (2016). Our transfer function had an observed versus predicted elevation r2 of 0.92 SWLI and a root mean squared error of prediction of 6.95 SWLI.

As with the foraminifera, we applied the modern analog technique to the fossil diatom assemblages using dissimilarity coefficients (minimum distance to closest analog, using the squared chord distance as the distance metric, MinDC, on Fig. S6; Kemp and Telford, 2015) to test the degree to which the assemblages in the modern samples provide analogs for the fossil assemblages. Samples with coefficients lower than the 20th percentile were defined as good analogs, and samples with coefficients larger than the 20th percentile were defined as poor analogs (e.g., Horton and Edwards, 2006; Kemp and Telford, 2015). Of our 136 fossil samples, 38 (28%) had MinDC values greater than the 20th percentile, including at least one sample adjacent to all contacts sampled except contacts B and C (Fig. S6).

The results of our diatom transfer function analysis are generally consistent with the results of the foraminiferal transfer function analysis in showing mostly minimal (<0.1 m) changes in RSL across contacts A through H (Fig. S6). The greatest change between means of elevation reconstructions for good analog samples across a contact suggests ∼0.12 ± 0.32 m of submergence across contact C. However, 30% of the mean elevations for diatom samples plot >0.1 m above foraminiferal reconstruction means for samples from the same levels in core S, whereas 34% plot >0.1 m below. Of course, reconstruction errors for both groups at the same levels overlap by a minimum of 20% (Figs. 9 and 10; Fig. S6).

We attribute differences between the elevations reconstructed with our diatom transfer function compared with those with the foraminiferal transfer function to the lack of a well-tested diatom transfer function for this region, comparable to those used in Alaska (Watcham et al., 2013; Shennan et al., 2016) or the Bayesian foraminiferal transfer function of Kemp et al. (2018). Our uncertainty in the accuracy of the diatom reconstructions is partly a reflection of the hundreds of diatom species, many with broad and(or) uncertain environmental preferences, that make up the diverse assemblages typical of Cascadia tidal sequences. Such assemblages have limited the degree to which modern diatom assemblages can be used as good analogs for fossil assemblages at Cascadia (e.g., Nelson et al., 2008). For example, at the Niawiakum River (a modern diatom site of Sawai et al., 2016) in the Willapa Bay region of southwest Washington, because most of Hong’s (2019) diatom samples from above and below six earthquake-subsided wetland contacts contained too many species to have good analogs in her modern data set, she grouped species with similar abundances across elevation to develop a diatom transfer function that resulted in improved analogs. Nevertheless, her earthquake-subsidence estimates for the six contacts using elevations reconstructed with the improved transfer function were significantly lower than other estimates of subsidence at the same and similar sites in the region. Because we are uncertain about the accuracy of our diatom transfer function reconstructions (Fig. S6), at the Siuslaw River, we relied on abundance data for the most common diatom taxa to qualitatively assess paleoenvironmental change across contacts in core S (Fig. 11; Table S4).

None of our stratigraphic evidence from beneath the marshes of the Siuslaw River estuary suggests meter-scale coastal subsidence during megathrust earthquakes, such as that reported for the greatest earthquakes at some sites (Atwater and Hemphill-Haley, 1997; Milker et al., 2016; Kemp et al., 2018). Do any of the 12 Cox Island contacts (Figs. 5 and 6) potentially record subsidence during a megathrust earthquake? To answer this question, we summarized lithologic contact characteristics and our correlation of contacts across the island and noted what foraminiferal faunas and diatom floras suggest about changing environments across contacts. Using fossil foraminifera, we then tested our inferences about contacts by reconstructing amounts of submergence (RSL rise) across 6 of the 12 contacts using the Bayesian foraminiferal transfer function of Kemp et al. (2018).

In assessing each contact, we relied on the most comprehensive review of criteria for identifying earthquakes in tidal wetland sequences: Shennan et al. (2016) expanded the original criteria of Nelson et al. (1996a) to include new methods, much greater detail in application of criteria, and additional criteria. Although they restated the importance of the original criteria of (1) the lateral extent of peat-mud contacts, (2) the suddenness of the change in environment across contacts, and (3) quantitative estimates of the amount of elevation change across contacts, Shennan et al. (2016) emphasized that such evidence should be consistent among multiple locations within the same coastal site. Early qualitative assessments, based primarily on inferences about where in the tidal zone particular lithologies typically form, concluded that changes in lithology suggesting coseismic elevation changes of roughly 0.5 m or less were difficult to distinguish from similar lithologic changes produced by nonseismic processes (e.g., Nelson, 1992b; Nelson and Kashima, 1993; Nelson et al., 1996a). More recently, Shennan et al. (2016) showed that, for Alaskan peat-mud contacts with reconstructed RSL rise or fall of <0.5 m, summary diagrams of diatom salinity preference have only a 50% chance of showing the correct trend of RSL between two samples across a contact. However, with sufficient quantitative, redundant, and consistent data for contacts at multiple locations at a site, Shennan et al. (2016) suggested that application of the above criteria may support an effective detection threshold for earthquake uplift or subsidence as low as 0.1–0.2 m.

By the criteria of Shennan et al. (2016), our evidence from the Siuslaw River estuary falls below the detection threshold for megathrust earthquakes because we lack quantitative microfossil evidence from more than a single core and nearby outcrop, and because no data suggest coseismic subsidence greater than 0.5 m (∼20% of the great diurnal tidal range near Cox Island). However, in our evaluation of each contact, we also considered the other criteria of Shennan et al. (2016), such as probable tsunami deposits above contacts and comparisons of Bayesian probability models for the times of earthquakes at the Siuslaw River with those for earthquake evidence at sites to the north and south.

Contact A

Contact A is more distinct and has been identified over a larger area than any other contact. Its distinctness is primarily the result of the sand, muddy sand, or sandy mud that caps muddy peat to peaty mud in many cores. In core S, a 2–3 cm bed of clean, very fine to fine sand, which grades upward into sandy silt, abruptly overlies slightly muddy peat at contact A, whereas in core K, indistinct 2- to 3-mm-thick laminae of silty sand above the contact fine upward into laminae of sandy silt, suggesting multiple depositional pulses. Three other gouge cores along transect A (Fig. 4) showed 2–4 cm of clean, very fine to fine sand overlying contact A. As contact A is one of only three contacts at Cox Island capped with sandy sediment (A, C, and Db; Figs. 5 and 6), and its sandy sediment thins upriver as well as landward, we infer that it is more likely a tsunami deposit than an extreme river flood deposit (e.g., Wells, 1947).

In the ∼40% of cores along transects A, B, and C that lack sandy sediment above contact A, peaty mud or mud commonly overlies muddy peat or peat, suggesting—at most—a few decimeters of RSL rise (e.g., Nelson et al., 1996a; Shennan et al., 2016). Where described in transect cores, the peaty sediment beneath contact A is typically 10–20 cm thick, with its upper contact much sharper than its lower contact. Contact A is mapped in cores 20–110 m apart for 200–400 m along the transects (Fig. 7). Although we lacked descriptions of cores over distances of as much as 400 m (Fig. 4), we are confident of our correlation of contact A from transect A (core 6, Fig. 5) to outcrop 2, a distance of 1000 m. In contrast, the interbedded peat, muddy peat, and peaty mud with largely gradual contacts in the upper 2–3 m of the reconnaissance gouge cores on two core transects in South Inlet (Figs. S2 and S3) and one transect near the mouth of the North Fork of the Siuslaw River (Fig. 2) showed a potential correlative of contact A in only 4 of the 16 cores described (one at 0.75 m depth in core 11; Fig. 3; Fig. S2).

To learn more about the change in environment marked by contact A, we analyzed assemblages in 15 foraminiferal samples (using the original data of Hawkes et al., 2011) and 42 diatom samples above and below the abrupt (1 mm) contact at 0.56 m in core S. The contact separates a muddy peat with 7.5YR hues from overlying sand, silty sand, and muddy peat (Figs. 5, 9A, and 9B). The upper four of the five foraminiferal samples below the contact are dominated by Balticammina pseudomacrescens (25%–55%), Trochammina inflata (12%–28%), Jadammina macrescens (10%–24%), and Haplophragmoides sp. (7%–43%), reflecting a middle to high marsh environment (Table S2; e.g., Hawkes et al., 2010; Engelhart et al., 2013a; Milker et al., 2015a). Hawkes et al. (2010) reported a similar assemblage above MHHW on the transect of modern foraminifera studied on Cox Island 600 m southwest of core S (transect M, Fig. 4). In the sand and silty sand above contact A, three samples were barren of foraminifera, but the lowest sample had a low-concentration assemblage dominated by B. pseudomacrescens (48%) and T. inflata (32%), perhaps as a result of mixing of foraminifera from the peat into the overlying sand (e.g., Milker et al., 2016). Three other low-concentration (33–39 tests; Table S2) assemblages from the muddy peat above the sand consisted almost entirely of J. macrescens (85%–100%), as did the three samples in the muddy peat above them (2–7 tests/mL).

Kemp et al. (2018) compared submergence (inferred to be subsidence) reconstructions for contact A using the fossil foraminiferal data of Hawkes et al. (2011) from core S as part of their development of a new Bayesian transfer function. Hawkes et al. (2011) had used these same data to reconstruct 0.4 ± 0.6 m (errors on all subsidence values at 2σ) of coseismic subsidence across contact A, but that reconstruction used an early non-Bayesian transfer function that was hampered by five no-modern-analog assemblages above the contact (no modern sample in the Hawkes et al. [2010] database was a good analog for the fossil assemblages). Using a much larger modern data set than Hawkes et al. (2011; 393 samples vs. 91 samples), Kemp et al. (2018) reconstructed subsidence for contact A with their non-Bayesian transfer function (0.1 ± 1.0 m), their Bayesian function with no lithologic priors (0.3 ± 0.8 m), and their Bayesian function with lithologic priors (0.5 ± 0.8 m). However, Kemp et al. (2018) assigned lithologic priors to the samples for contact A following the simplified lithology for core S shown in Hawkes et al. (2011, their fig. 3d). Our more detailed lithologic description of the same core (Fig. 5) suggests that four of the samples above the barren samples above contact A probably formed between MHW and MHHW and, therefore, might be placed in either of the two lithologic prior groups of Kemp et al. (2018; Table S2). Our reanalysis using the new Bayesian function with lithologic priors (with minor corrections to some of the fossil data) gives 0.5 ± 0.8 m of subsidence with the lithologic priors used by Kemp et al. (2018), but 0.1 ± 1.0 m if the four samples are placed in Kemp et al.’s (2018) higher lithologic prior group (Figs. 9A and 9B; Table S2). Thus, considering the uncertainty in selecting the most appropriate priors for samples of muddy peat that were probably deposited between MHW and MHHW, the range in subsidence obtained with different transfer functions, and the low concentration of foraminifera in all samples above contact A, we conclude only that—based on the foraminifera—subsidence across contact A was probably <0.5 m.

The modest increases in the abundance of freshwater diatom taxa and decreases in brackish taxa across contact A are more consistent with a slight decrease in salinity rather than a significant increase (Figs. 9 and 11; Tables S3 and S4). Species with a low-salinity preference, such as Cosmioneis pusilla (5%–15%), Pinnularia lagerstedtii (7%–12%), and Luticola mutica (8%–22%), and brackish species, such as Navicula cinta (6%–12%), dominate assemblages of the muddy peat beneath the contact. The former species are some of the most common diatoms in the high marsh along the modern transect on Cox Island and other Oregon transects studied by Sawai et al. (2016) (transect M on Fig. 4). Immediately above the contact, these species are almost absent in the sand, which is dominated by Planothidium lanceolatum, a species common in flowing water, Rhoicosphenia abbreviata, a freshwater-brackish species, and Cocconeis placentula, an epiphytic freshwater-brackish species (Fig. 11). In the peaty mud 5–10 cm above the contact, the above low-salinity-preference species again become dominant, along with greater percentages of freshwater species and a freshwater-brackish species, Caloneis bacillum (12%–30%).

Because subsidence across contact A reconstructed with the foraminiferal transfer functions is variable, and the diatom assemblages across the contact are most consistent with minimal changes in tide levels (Fig. 11; Fig. S6), subsidence across the contact was below the 0.5 m detection threshold for distinguishing coseismic subsidence from nonseismic origins, especially without additional reconstructions across correlative contacts elsewhere in the estuary (Shennan et al., 2016).

If there was minimal subsidence across contact A, could its sandy capping bed be a river flood deposit rather than a tsunami deposit? Diatom assemblage changes are inconsistent with a marine origin (Fig. 11). However, based on (1) the lithologic distinctness and wide distribution of contact A on Cox Island, (2) its shallow depth (∼0.5–0.6 m), (3) changes in foraminiferal assemblages and lithology across it that suggest at least 0.1 m of subsidence, and (4) laminae of clean sand suggestive of tsunami pulses within its capping sandy bed, we favor a tsunami origin for the bed. For this reason, we correlate the sandy bed and its lower contact (A) with evidence for the 1700 CE (250 cal yr B.P.) earthquake and its tsunami along much of the subduction zone (e.g., Nelson et al., 1995, 2006; Witter et al., 2003; Atwater et al., 2004, 2005; Graehl et al., 2014; Valentine et al., 2012; Milker et al., 2016; Hutchinson and Clague, 2017; Hong, 2019; Padgett, 2019). Five 14C ages from core S and outcrop 1 (Fig. 6; Table 2) are similar to ages for evidence of the 1700 CE earthquake at many tidal sites from Washington to northern California. Although the five ages are consistent with our correlation, such young 14C ages on detrital macrofossils only show that contact A is less than ∼500 yr old (e.g., Kemp et al., 2013). Geophysical models of the rupture during the great earthquake of 1700 CE (discussed below) are consistent with our correlation in that they show coseismic subsidence of ∼0.2 m near the Siuslaw River estuary (Wang et al., 2013, their fig. 8; Wirth and Frankel, 2019, their fig. 3c).

Contact B

Because contact B is much less distinct, more difficult to correlate along core transects, and much less widespread than contact A (Fig. 7), we are uncertain of its origin. Along transect A (Fig. 4), two contacts with peaty mud overlying peat or muddy peat occur 15–50 cm below contact A in some cores, but the single such contact in this stratigraphic interval in other cores makes correlations uncertain. Where both contacts are present, the upper contact (B) is always sharper (9 ± 13 mm, 1σ error), although little sharper than the lower contact of the peaty unit beneath contact B (Fig. 7). Although we show its possible correlation on Figure 5, in most cores along transects B and C, contact B is either absent or indistinguishable among several gradational contacts between muddy peat and peaty mud units in the upper 1.5 m interval of the cores.

Diatom assemblages in samples above and below contact B in core S suggest minimal changes in salinity that might reflect changes in tidal environments across the contact (Fig. 11; Tables S3 and S4). Counts of two of the three dominant species decrease slightly across the contact (C. pusilla, upward change from 1%–33% below to 0%–10% above; P. lagerstedtii, 6%–9% to 7%–10%; L. mutica, 6%–27% to 5%–21%), whereas percentages of the most common freshwater-brackish species change little (Navicula cryptotenella, 1%–7% to 1%–6%; N. cincta, 1%–9% to 5%–12%). As the correlation and origin of contact B were uncertain, we did not sample it for foraminifera.

An oldest minimum age of 210 ± 25 14C yr B.P. for contact B in our final age model gives an interval of 354–165 cal yr B.P., which overlaps considerably with our interval for contact A (Fig. 8; Table 2). Even if older by a century or two than contact A, contact B’s age interval does not overlap with the intervals for any published evidence of pre–1700 CE earthquakes or tsunamis in Oregon (Fig. 12). It may record a small, gradual, very localized change in tide levels resulting from changes in river-channel or estuary configuration combined with gradual sea-level rise rather than sudden coastal subsidence during a great earthquake (e.g., Nelson et al., 1996b).

Contact C

Contact C is too indistinct to map with certainty along all of transect A, but it is more distinct and continuous (>400 m) along transects B and C (Fig. 7). Contact C is particularly distinct where sandy mud or muddy sand cap peat or muddy peat (in 47% of cores; e.g., Fig. 5), where 2 cm of slightly silty very fine sand overlie silty peat in three gouge cores along transect A (Fig.3), and near the west end of outcrop 1, where a 1.2-m-long spruce root is rooted in the peat below it (Fig. S5). Elsewhere, contact C separates peat or muddy peat from overlying mud or slightly peaty mud with mean contact thicknesses of 4–9 mm (Fig. 7). Although the sandy beds above the contact are widespread enough—even 100–200 m inland from the river—for us to infer deposition by a tsunami, the muddy sand without distinct laminae suggestive of tsunami pulses may also have been deposited by a river flood.

The minimal changes in foraminiferal assemblages across contact C at outcrop 1 seem inconsistent with the lithologic change from a slightly muddy peat overlain by a silty sand or sandy mud grading upward into a slightly organic-rich mud and even more inconsistent with our Bayesian transfer function reconstruction of submergence across the contact of 0.4 ± 0.6 m (Fig. 9D; Table S2). High-abundance foraminiferal assemblages above and below the contact (depth of 91 cm; Fig. 9D; Table S2) are dominated by B. pseudomacrescens (40%–61%), J. macrescens (18%–35%), and Haplophragmoides sp. (12%–24%), with the only significant change in assemblages across the contact being low numbers of Trochamminita sp. (8%–14%) below it and T. inflata (2%–4%) above it (Table S2). Such assemblages are typical of Oregon middle marshes above MHW (Hawkes et al., 2010), but if the mud above the contact were deposited in a low marsh below MHW, we would expect the mud to host significant percentages of Milliammina fusca, by far the most common foraminiferal species of low marsh and tidal flat environments in the region (Kemp et al., 2018).

Although the lithologic contrasts across contact C in core S are less distinct than those at outcrop 1—where a slightly muddy peat is overlain by a very muddy peat with coarse silt near its base—diatom assemblages in samples across contact C in core S (depth of 95 cm; Tables S3 and S4) are consistent with a salinity increase that could reflect decimeters of submergence (Fig. 11; Fig. S6). The abundances of some species, such as Navicula tenelloides (9%), are largely unchanged across the contact, whereas L. mutica (upward change from 12%–24% below to 1%–11% above), Navicula rhynchocephala (0% to 6%–13%), N. cryptotenella (1%–4% to 4%–8%), and especially Navicula gregaria (4% to 43%–59%), significantly increase. Low numbers of the marine species Rhaphoneis surirella (3% to 11%) below the contact show a fourfold upward increase across the contact as well.

Contact C is unique in being younger than other pre–1700 CE tidal wetland stratigraphic contacts inferred to record subsidence during megathrust earthquakes in coastal Oregon and Washington (e.g., Atwater and Griggs, 2012, p. 22; Garrison-Laney, 2017). An age (310 ± 30 14C yr B.P.) on seeds of an herb (cf. Atriplex sp., commonly found in the upper high marsh) from the peat below contact C at outcrop 1 closely overlaps the age (335 ± 30 14C yr B.P.) on rings 7–8 of the spruce root rooted in the peat (Fig. 6; Table 2), suggesting these ages, along with another wood root age (310 ± 15 14C yr B.P.), are close maximum ages for contact C. Using the youngest of the three maximum ages in our age model gives an interval of 463–283 cal yr B.P. for contact C. From the concordance of the three ages on high-quality samples, we infer that contact C is a century or two older than contact A. Darienzo et al. (1994) reported bulk 14C ages in this age range on muddy peat and peaty mud in cores from seven estuaries in northern Oregon, but such ages could be hundreds of years older or younger than the sampled levels in the cores (e.g., Nelson, 1992a). The only widely reported evidence attributed to a megathrust earthquake that overlaps significantly with this time period is marine turbidite T2 of Goldfinger et al. (2012) (Nelson et al., 2006; Hutchinson and Clague, 2017; Fig. 12 here).

Evidence is insufficient to infer coseismic subsidence for contact C followed by deposition by a tsunami. The change in environment indicated by diatoms and the submergence reconstructed with the foraminifera are probably below the 0.5 m detection threshold for coseismic subsidence. In any case, the reconstruction is based on a single sample; other samples showed no change in tide level across the contact (Fig. 9D). Based on its >400 m extent along transects B and C, its upward lithologic change from a muddy peat to an organic-rich mud, its overlying muddy sand, and the changes in diatom assemblages consistent with an increase in salinity across it, contact C might coincide with regional subsidence during a megathrust earthquake. However, the lack of evidence for an earthquake about this time at other coastal sites where evidence for earlier earthquakes is quite distinct makes a coseismic origin for contact C unlikely. If its overlying sandy bed were deposited by an extreme river flood, it may record a small, gradual, very localized change in tide levels resulting from changes in river-channel or estuary configuration combined with gradual sea-level rise.

Contacts Da and Db

Along transect A, a single contact D sharply separates muddy peat from overlying rooted mud, but along much of transects B and C, the underlying peaty unit consists of two peaty beds separated by a 5- to 15-mm-thick, light-gray, rooted, silt lamina. We labeled the contact below the lamina contact Db and the much more gradational contact above it contact Da (Figs. 5 and 6). Our correlation of the 400 m extent of the two contacts at the stratigraphic position of contact D along transects B and C relies primarily on the most distinct contact, Db (Fig. 7). The contact correlated in most cores along transect A is probably also Db, but the absence of the silt lamina in this area makes correlations less certain. Three of the 21 cores in which we identified contact Db showed sandy rooted silt or sandy peaty mud above the contact. In core S, three laminae overlie the well-humified high marsh peat with 7.5YR color hues beneath contact Db: 7–9 mm of slightly sandy peaty silt, 6–10 mm of silty fine sand, and 8–10 mm of silt with coarse fragments of organic debris in the upper 4 mm (Fig. 5). Such distinct laminae, each of which fines upward, are more typical of pulses of deposition during tsunami inundation than of river flood deposits. As for contact A, we infer that the sandy laminae above contact Db were more likely deposited by a tsunami than by a river flood. If so, the extensive silt lamina above contact Db in many cores may record suspension deposition from a tsunami surge between inundation and return flow.

Low concentrations (1–17 tests/mL; Table S2) of foraminifera in the two samples above contact Db in core S are consistent with a tsunami origin, whereas assemblages in two samples below the contact and three samples above the low-concentration samples suggest middle to high marsh environments returned soon after deposition of the sandy beds (Figs. 5 and 9C). The two samples from the peat below the contact contain B. pseudomacrescens (36%–57%), T. inflata (22%–29%), Haplophragmoides sp. (11%–16%), J. macrescens (1%–10%), and M. fusca (6%–9%) (Table S2), similar to the assemblages near MHHW on the modern transect of Hawkes et al. (2010) (transect M, Fig. 4). Except for the sample with 20% M. fusca just above the sandy beds, the other two samples above the contact have percentages similar to those for the same species below the contact. Although we used the sample with 20% M. fusca above the sandy beds to reconstruct the RSL change across the contact, the Bayesian transfer function suggests only a small RSL rise across it: 0.1 ± 0.6 m (Fig. 9C). As contact Da was less distinct and much more difficult to correlate than contact Db, we did not sample it for foraminifera.

Diatom assemblages in the three diatom samples above and three samples below contact Da suggest, perhaps, a slight increase in salinity (Fig. 11). The abundant species C. pusilla (33%–5%) and L. mutica (20%–7%) decreased significantly in the sample 1 cm above the contact, but then increased to similar levels in the two higher samples. In contrast, Nitzschia sigma (1%–8%) and Tabularia fasciculata (9%–12%) increased in the sample above the contact and then decreased upward, as did Gyrosigma eximium (1%–9%). The dramatic increase in the brackish to marine species Paralia sulcata (0%–68%) above the contact may indicate a temporary influx of more saline water, perhaps during a storm surge or tsunami.

For contact Db, L. mutica, C. pusilla, and P. lagerstedtii retained their dominance in the four samples above and five samples below the contact, although their abundance dropped considerably in the samples from the probable tsunami deposit 2–4 cm above the contact (Fig. 11). Rhoicosphenia abbreviata (upward change from 0% below to 5%–17% above) makes a sudden appearance above the contact, as do the freshwater species P. lanceolatum (0% to 10%) and Gomphonema parvulum (0% to 4%). However, the large upward increases in the percentage of the genus Thalassiosira sp. (6% to 14%–19%) across the contact in these three samples suggests an influx of plankton. The lack of other brackish species in the samples above those containing the plankton suggests the influx was due to a short-lived storm surge or tsunami rather than to decimeters of coseismic subsidence. This inference is consistent with the foraminiferal transfer function reconstruction of 0.1 ± 0.3 m of submergence across contact Db (Fig. 9C).

Although the ages of contacts Da and Db are each constrained with minimum as well as maximum ages, the minimum ages provided by herb roots are probably much younger than the times the contacts formed (Table 2). The youngest seeds sieved from peat below the contacts probably give the closest estimates of contact age. Although the two contacts are only 0.2 m apart, their age intervals barely overlap because of the two-and-a-half century difference in their minimum ages: Da, 676–511 cal yr B.P. and Db, 790–670 cal yr B.P. Three other ages from the peat below each contact are older, suggesting all are on detrital materials.

The age interval for contact Da overlaps with only the youngest portions of some modeled age distributions for earthquake and tsunami evidence of about this age at other central Cascadia sites (Fig. 12). Diatom assemblages show little evidence for significant or long-lasting environmental change across contact Da, and our difficulties in mapping it beyond core S and outcrop 1 prevent determination of its origin. It probably records a small, perhaps very localized change in tide levels resulting from changes in river-channel or estuary configuration, perhaps during a storm surge or tsunami.

Although contact Db is considerably more distinct and extensive (400 m) than contact Da, neither diatoms nor foraminifera suggest a significant, long-lasting change in environment or RSL across it. Its age distribution, however, overlaps considerably with distributions for evidence at six other coastal sites on Figure 12, as well as turbidite T3 of Goldfinger et al. (2012). If our inference of tsunami deposition for its overlying beds is correct, it records a tsunami from a megathrust earthquake that produced little measurable subsidence at Cox Island.

Contacts Ea through H

Contacts Ea, Eb, Fa, Fb, G, and H are the most distinct of a series of contacts separating 1- to 4-cm-thick beds of peat, muddy peat, peaty mud, and mud in the lower 1.5 m of (compacted) core S (∼2.0–4.2 m depths on the uncompacted version of core S in Fig. 5). As for contact D (Da, Db), we correlated contact E as a single contact along transect A and parts of transects B and C. However, we subdivided it into contacts Ea and Eb in cores near the intersection of transects B and C. Only the lower contact, Eb, was identified in outcrop 1 (Fig. 6). Likewise, two contacts (Fa and Fb) of peaty mud over peat or muddy peat were distinct in core S, but only contact Fb was described from the outcrop. However, correlations of contacts below contact D (Da, Db) were less certain than for younger contacts because below this depth, most cores display a succession of alternating muddy peat, peaty mud, and peat units of varying thickness with only a few, thin intervening sections of rooted mud (Fig. 5). Some of the 1- to 2-cm-thick beds of peat may be detrital. Although we are confident in the identity of contacts in core S (Fig. 5), our detailed core descriptions record a variable series of muddy peat and peaty mud laminae and beds above and below each contact. Only 10–12 cores reached these deeper contacts, and contacts Ea, Fa, and H were not identified at the outcrop.

For the peaty units below contact D, some upper contacts were sharper than lower contacts, but for other contacts, thicknesses were similar, and the range in thickness for all contacts was great (1σ for thicknesses ranged from 50% to 100%; Fig. 7). In most cores, the deeper contacts lacked consistent distinctive characteristics, such as erosional topography, strong contrasts in lithology, or caps of coarse silt or sandy mud that could be used to distinguish them from similar upper or lower contacts. In both the adjacent vibracores at site S, contacts Ea, Eb, and G were marked by middle to high marsh peat with 7.5YR color hues sharply overlain by mud or peaty mud, whereas contacts Fa, Fb, and H were less distinct with muddy middle marsh peat overlain by peaty mud. Our correlations of contacts E (Ea, Eb), F (Fa, Fb), G, H, and I relied on relative depth in the sequence and the sharpest contacts at the tops of the most peat-rich units, but other correlations of contacts are possible. Radiocarbon ages supported our correlations of contacts Eb, Fb, G, H, and I between core S and outcrop 1 (Fig. 6). Contacts E and F on transect B, and contacts F, G, and H on transect C were traced for >300 m, but we correlated most other contacts <200 m (Fig. 7). Even the greatest lithologic contrasts across the deeper contacts—peat overlain by mud or rooted mud—were typical of <1 m of rapid RSL rise (e.g., Shennan et al., 2016; Horton et al., 2018), and most lithologic contrasts suggested less submergence.

The succession of interbedded peat, muddy peat, and peaty mud with largely gradual contacts in the upper 2–4.4 m of the reconnaissance gouge cores described in South Inlet and the North Fork of the Siuslaw River showed only indistinct lithologic contrasts, typical of less than a few decimeters of RSL change (Fig. 3; Figs. S2, S3, and S4). Except in cores along the north edge of South Inlet and in the upper parts of a third of the cores from the North Fork of the Siuslaw River, contacts were less distinct and continuous than on Cox Island. Without many more 14C ages from cores at these sites, correlations to contacts on Cox Island are speculative.

Foraminiferal analyses of samples across contacts Ea, Eb, Fa, Fb, G, and H showed diverse results reflecting environments of variable salinity and exposure. Although sampled lithologies were similar to younger units characterized by foraminiferal assemblages representative of tidal environments, the absence of foraminifera in many samples near these contacts suggested extended periods of largely freshwater deposition. Three samples across contact Ea in core Sb were barren, as were four samples in core Sb and four samples from outcrop 1 across contact Eb (Table S2). In contrast, only one sample among eight across contact Fa in core Sa was barren (2 had <14 foraminifera). Nevertheless, changes in principal species typical of the high marsh in the four samples spanning contact Fa were minimal: B. pseudomacrescens (29%–50% below to 18%–50% above), Haplophragmoides sp. (35%–50% below to 35%–41% above), and J. macrescens (5%–21% below to 10%–34% above). Changes in assemblages in the two nonbarren samples of the four samples across contact Fb at outcrop 1 were even less pronounced, with the only difference above and below the contact being a change from 20% to 8% T. inflata (Table S2). Of the 13 samples from core Sb and four samples from the outcrop across contact G, 10 were barren, and another three had low concentrations of foraminifera. Although two of the other four samples spanned contact G, both consisted of 100% Milliammina petilla (a tidal-flat species), and so they were not useful for reconstructing RSL change. Similarly, of the 11 samples spanning contact H in core Sb, all samples near the contact were barren or low-concentration samples (Table S2).

Reconstructions of RSL rise with the Bayesian transfer function across the only two of these contacts with apparently in situ foraminiferal assemblages immediately above and below the contacts gave 0.2 ± 0.3 m for contacts Fa and Fb (Figs. 10A and 10B; Table S2). Because of the few samples with >30 foraminifera and the similarity of the assemblages on either side of the contacts (for example, neither of the samples above the contacts contained M. fusca), it was difficult to assess the accuracy of these reconstructions. A change from middle to high marsh peat to peaty mud across the contacts is consistent with a few decimeters of submergence, but the reconstructed submergence for both contacts is well below the 0.5 m threshold of detection for coseismic subsidence.

Unlike samples near higher contacts in the core, samples above and below contacts Ea and Eb contained fewer diatom species with a low-salinity preference and very few brackish species (Fig. 11). The main exceptions were Gomphonema gracile, which increased significantly above both contacts (Ea, upward change from 4%–6% below to 2%–29% above; Eb, 0%–2% to 8%–12%), as do R. abbreviata (9%–11% to 19%–24%; Eb 9%–11% to 10%–22%), C. placentula (Ea, 1%–2% to 8%–18%; Eb, 0% to 3%–4%), and Rhoicosphenia linearis, which increased above contact Eb (1% to 7%–14%). The freshwater species Encyonema minutum (2%–16%), Eunotia praerupta, G. parvulum, and P. lanceolatum all increased significantly above contact Ea. Across contact Eb, P. lanceolatum remained unchanged, whereas Gomphonema gracile increases (2% to 8%–21%) and G. parvulum decreases (15%–26% to 0%–18%). Relatively few valves of brackish and marine species were counted in samples near contacts Ea and Eb (Fig. 11). Although these assemblage changes are modest, they are consistent with environments of primarily low salinities becoming fresher across contact Ea. The absence of foraminifera in all samples near contacts Ea and Eb is consistent with the primarily freshwater environments indicated by the diatoms. Thus, neither foraminifera nor diatom assemblages near contacts changed sufficiently to infer a coseismic origin for contacts Ea and Eb.

Much like contacts Da and Db, the most abundant species above and below contact Fa in core S are L. mutica (upward change from 0% below to 30% above), P. lagerstedtii (1% to 28%), and C. pusilla (1% to 47%), species common in the high marsh (Sawai et al., 2016). Although the valves of these species were partly replaced by those of the river species P. lanceolatum in the sample directly above the contact (0% below to 20% above), in higher samples the same high marsh taxa return in significant numbers. As for contact Db, the six-fold increase (0% below to 63% above) of the planktonic genus Thalassiosira spp. in the peaty mud above the contact is consistent with a short-lived influx of marine water, perhaps during a storm surge, before species common in the high marsh returned. In core S, contact Fb was less distinct and more gradual (a 4-mm-thick contact) than most other contacts and so was not sampled for diatoms.

Diatom assemblages above and below contact G were dominated by species with a low salinity preference and, to a lesser extent, by a few freshwater species. Species that increased significantly across the contact include Navicula pseudolanceolata (0% below to 8% above) and G. gracile (3% to 17%), whereas C. pusilla (24% to 1%) and P. lagerstedtii (26% to 0%) decreased. P. lanceolatum, common in flowing water (Patrick and Reimer, 1966), increased significantly above the contact (0% to 14%), and the additional freshwater species Gomphoneis mammilla (0% to 14%) and Encyonema minutum (0% to 9%) made an appearance above the contact. However, other freshwater species, such as Ulnaria ulna, did not change across the contact. We infer that such assemblages may reflect a fluctuating tidal environment more strongly influenced by river flooding than changes across higher contacts, but without any significant, long-lasting changes in environment. The 13 barren and low-concentration foraminiferal samples (eight containing freshwater thecamebians) from contact G to 15 cm below it in core Sb (Table S2) are consistent with a fluctuating but primarily freshwater to slightly brackish environment. At least once prior to contact G, however, a high marsh was regularly inundated by brackish water at the site: 10 cm below contact G in core Sb, two samples contained a typical high marsh foraminiferal assemblage. The two other samples with adequate concentrations of foraminifera that spanned contact G in the same core contained only 42–93 tests of the tidal-flat species Milliammina petula, which is inconsistent with the diatom assemblages near the same contact in adjacent core Sa. Again, neither foraminifera nor diatoms suggested sudden submergence on a scale typical of coseismic subsidence identified in other Cascadia coastal sequences.

In contrast to contact G, diatoms across contact H showed an apparent change from a primarily freshwater environment to a brackish environment (Fig. 11). This was reflected by species with a preference for freshwater (e.g., Stauroneis phoenicenteron [5%], P. lanceolatum [4%], G. parvulum [8%], Eunotia praerupta [16%], and E. minutum [2%]) in samples below the contact changing to mostly freshwater-brackish species (e.g., N. sigma [11%], N. cryptotenella [5%], N. rhynchocephala [5%–8%], Navicula libonensis [1%–9%], and P. lagerstedtii [1%–8%]) above the contact. The brackish species Caloneis westii also increased above the contact (0%–29%). The eight barren and low-concentration foraminiferal samples (five containing freshwater thecamebians) from near contact H in core Sb (Table S2) are consistent with a fluctuating but primarily fresh to slightly brackish environment. As with the samples near contact G, however, foraminiferal assemblages typical of Oregon high marshes in one sample below and two samples above contact H (Table S2) show that brackish environments with regular tidal inundation occasionally characterized the site of core S. Although diatom data cannot preclude increasing salinity as a result of a few decimeters of coseismic subsidence, the data are well below the threshold of detection for such an event (Shennan et al., 2016).

OxCal age-model intervals for contacts E through H relied on our assessment of the most accurate maximum and minimum ages from above and below each contact. Two to five older maximum ages on detrital materials for each of these contacts are consistent with the modeled age intervals (Table 2). Concordant maximum ages on two high-quality detrital samples from the peat below contact Ea suggest that, at 940–770 cal yr B.P., it is less than a century older than contact Db. Similarly, the youngest of three concordant maximum ages from the peat below contact Eb yielded an interval of 1207–922 cal yr B.P. for that contact. Assuming the much younger age on seeds of 1330 ± 30 14C yr B.P. is an outlier, the youngest maximum and oldest minimum ages for contact Fa gave an interval of 1455–1308 cal yr B.P. The youngest age on a stem base and rhizome below contact Fb limited its age to 1512–1398 cal yr B.P. If this rhizome and another at about the same depth grew down into the underlying peat from above the contact, abundant herb seeds in the peat would still limit the contact’s age to 1535–1423 cal yr B.P. Unless an herb rhizome in growth position 5 mm above the contact is detrital, the youngest maximum and oldest minimum ages for contact G yielded an interval of 1636–1508 cal yr B.P. Ages on the oldest of two rhizomes and herb seeds from the peat below contact H show it to date from 1778 to 1588 cal yr B.P. If the rhizomes do not provide minimum ages for the contact but instead are detrital, contact H would be only slightly older than contact G (1674–1557 cal yr B.P.).

Although the age distributions for contacts Ea, Eb, Fa, Fb, G, and H partially overlap distributions for earthquake and tsunami evidence at other sites on Figure 12, lithologic and microfossil evidence is inconclusive about possible coseismic origins for any contacts. The probability distribution for contact G is a good match for the earthquake ca. 1550 cal yr B.P., which has much evidence indicating >0.5 m of subsidence farther north (Fig. 12). Although contact G’s diatom assemblages may be consistent with a few decimeters of sudden submergence, the diatom data are well below the threshold of detection for distinguishing local river or tidal changes from earthquake subsidence (Shennan et al., 2016).

Contact I

Although its extent is limited and uncertain because it was reached in only two to five cores along each transect (Fig. 7), the lithologic contrast across contact I was greater than for younger contacts below contact C. For example, in core S, contact I sharply (1–2 mm) separates a high-marsh peat with 5YR color hues from the overlying 6 cm interval of mud to peaty mud. However, the lower half of the mud is unusual in having a much higher proportion (close to half) of clay than does most estuarine tidal-flat mud. Such organic-rich fine mud is more typical of quiet-water lagoons than of tidal mudflats. The trace of very fine sand in the lower half of the mud may be evidence for river or tidal flooding soon after the contact formed. The 2.5-m-long cluster of spruce roots at this level at outcrop 1 shows that large trees were growing along the river about the time that contact I formed (Fig. 6).

Although the two foraminiferal samples immediately above and below the contact were barren, the one underlying sample and two of the four overlying samples suggested a large RSL rise across contact I (Table S2). The sample below consisted of 57% Haplophragmoides sp. with 43% B. pseudomacrescens, species typical of the high marsh; however, other expected high marsh species, such as T. inflata and J. macrescens, were absent. The two overlying samples contained 96%–100% A. mariae, a tidal-flat species that commonly occurs with M. fusca. The next higher sample hosted 47% M. fusca, a species commonly dominant in the low marsh, and 53% B. pseudomacrescens, a species characteristic of the high marsh (Hawkes et al., 2010; Milker et al., 2016). Other species that typically would occur with these species were absent.

Using the elevations reconstructed with the foraminiferal samples above and below contact I, the Bayesian transfer function gave 1.6 ± 0.8 m of submergence across the contact (Table S2), a rise 15% greater than any of the rises calculated by Kemp et al. (2018) for the 1700 CE earthquake at 15 sites spanning >400 km of the subduction zone. However, because none of the samples used by Kemp et al. (2018) to develop their transfer function contained high percentages of A. mariae, the two samples above contact I have no modern analogs in the Kemp et al. (2018) data set, making the reconstruction problematic. As A. mariae is a taphonomically resistant species (Goldstein and Watkins, 1999), perhaps other foraminiferal species originally deposited with it have decayed.

Although some changes in diatom assemblages occur across contact I, the dominant species on either side of the contact indicate fresh and low-salinity preferences. Below the contact, freshwater species, such as Pinnularia notabilis (9%), G. parvulum (24%), G. mammilla (8%), and E. praerupta (24%) were found to be dominant, although G. gracile (6%) and the aerophilic C. pusilla (7%) also occurred. Above the contact, several of the freshwater species decreased (although the abundance of G. mammilla doubled), and R. abbreviata (8%) and G. gracile (24%) appeared or increased. The brackish-marine species R. linearis also increased above the contact (1%–8%), and a brackish species, N. cincta, made an appearance (2%). All assemblages are inconsistent with the strongly brackish environment of a tidal flat suggested by the A. mariae fauna of the foraminiferal samples above contact I.

For this reason, and because the foraminiferal species assemblages in all but the highest sample, 8 cm above contact I, are not typical of Oregon tidal environments, we did not consider the foraminiferal submergence reconstruction for contact I reliable. Because the foraminiferal assemblages above and below contact I may not be representative of their environments of deposition, and the diatoms are consistent with an environment of low salinity changing to, at most, a slightly more brackish environment, we could not determine how much submergence, if any, is represented by the lithologic change across contact I. The high clay content with a trace of very fine sand and fresh and fresh-to-brackish diatoms in the lower 2 cm of the mud above the contact are consistent with flooding during the breaching of a lagoon. Although such breaching could occur during river flooding, transport of A. mariae tests by a storm surge or tsunami into a lagoon would also explain their presence in the mud above contact I. As A. mariae can withstand limiting environments with minimal light, oxygen, and salinity (Duijnstee et al., 2003), perhaps this species could even reproduce in a slightly brackish lagoon.

Our closest minimum and maximum ages for contact I place it in the interval 2128–1900 cal yr B.P. This assumes that the spruce stump rooted near contact I in the outcrop died about the time the contact formed and that a younger age on the outermost wood rings on one of its roots includes carbon that penetrated the root after death (Table 2). Four other ages on various materials from near contact I are consistent with this age interval. Studied stratigraphic sections at many sites do not sample events of this age, and only two other sites on Figure 12 have broad age distributions that overlap substantially with our distribution for contact I. Although lithologic changes across contact I are distinct, and spruce trees were rooted in a wetland soil horizon below the contact, diatom assemblages that conflict with the foraminiferal reconstruction of subsidence and the absence of age distributions for correlative earthquakes at most other sites make it highly uncertain whether or not contact I records earthquake subsidence.

None of the 12 studied contacts on northern Cox Island yielded criteria sufficient to show >0.5 m of coseismic subsidence and, therefore, a megathrust earthquake origin. Foraminiferal transfer function reconstructions across contact A in core S, which range from 0.1 ± 1.0 m to 0.5 ± 0.8 m (2σ errors), make the amount of subsidence across the contact uncertain (Figs. 9A and 9B; Table S2). However, based on its depth, wide stratigraphic extent, and probable overlying tsunami deposit, contact A probably correlates with more distinct evidence to the north and south dating from 1700 CE. In their preferred model of the complex rupture processes during great earthquakes along the subduction zone, Wang et al. (2013) inferred four patches of high-moment strain release along the 1700 CE rupture separated by areas of low-moment release, one within 20 km of the Siuslaw River estuary. However, as with the preferred model of Wang et al. (2013, their fig. 8), subsidence at the Siuslaw River estuary in the most consistent model presented by Wirth and Frankel (2019, their fig. 3c) is ∼0.2 m.

In a series of alternative models for the 1700 CE rupture using a more detailed three-dimensional, seismic-wave velocity model of the subduction zone, Wirth and Frankel (2019) distinguished between (1) large rupture patches of more uniform, low-to-moderate slip in the downdip region of the megathrust, including small, higher-slip subevents (earthquakes), which would have generated high-frequency energy (that strongly impacts buildings and infrastructure), and (2) broad, high-slip rupture patches that would generate low-frequency energy (most important for generating tsunamis). Within their series of magnitude 9.2 earthquake scenarios, Wirth and Frankel (2019) explained the apparent along-strike slip heterogeneity indicated by the large differences in coseismic subsidence at sites only tens of kilometers apart in central Oregon (Kemp et al., 2018) by inserting a high-frequency subevent (a magnitude 8.2 earthquake in the downdip part of the rupture) 70–150 km north of the Siuslaw River.

None of the evidence for other contacts at Cox Island, or at other reconnaissance gouge core transects in South Inlet and the North Fork of the Siuslaw River, is sufficient to conclude that they record subsidence during a megathrust earthquake. More detailed diatom or foraminiferal analyses at multiple sites (e.g., Shennan et al., 2016; Padgett, 2019) in the Siuslaw River estuary might show that some contacts record a few decimeters of rapid subsidence that might be attributed to earthquake subsidence. For example, contact Db has good evidence for a tsunami, and contacts G and H have age intervals that overlap with the intervals for one of the greatest earthquakes of the past 2000 yr at sites farther north. However, the predominance of diatom species with fresh or low-salinity preferences throughout core S, and its high proportion of barren and low-concentration foraminiferal samples, may reflect more local rather than regional RSL changes of <0.5 m during the past 2000 yr, as might be expected along a river with a Coast Range drainage basin of 2000 km2 that is 6 km upriver from the sea. The many peaty units with gradual contacts in reconnaissance gouge cores along transects in South Inlet and the North Fork of the Siuslaw River are consistent with this inference (Nelson, 1992b) (Figs. S2, S3, and S4). Alternatively, our lack of evidence for subsidence near Cox Island does not suggest that megathrust ruptures did not extend along the central Oregon coast during this period, only that they likely produced <0.5 m of coseismic subsidence.

Such conclusions, however, set useful limits on models of successive megathrust ruptures by suggesting that coseismic subsidence at this site was <0.5 m for megathrust earthquakes that may have ruptured this part of the subduction zone during the past 2000 yr. As all rupture models of slip heterogeneity along the subduction zone rely on microfossil-based reconstructions of coseismic land-level change, a <0.5 m limitation is valuable in assessing the region’s seismic hazard. For example, it could help limit the landward extent of megathrust rupture or suggest whether or not the locations of high-frequency subevents remain stable over multiple earthquake cycles (Wang et al., 2013; Gao et al., 2018; Wirth and Frankel, 2019). Such a limitation is also consistent with a long-lived structural feature influencing the distribution and strength of slip patches on megathrust ruptures along this part of the subduction zone, such as an offshore basin margin high or the seamounts being subducted northwest of the Siuslaw River (Fig. 1; Wells et al., 2003; Tréhu et al., 2012; Wang et al., 2013; Wang and Tréhu, 2016; Wirth and Frankel, 2019).

Lithostratigraphy and biostratigraphy beneath the marshes of Cox Island largely confirm the interpretations of Nelson (1992b), based on reconnaissance gouge coring, about the character of late Holocene relative sea-level (RSL) rise at the Siuslaw River estuary. More importantly, they limit geophysical models of Cascadia megathrust rupture during successive earthquakes by ruling out substantial coseismic subsidence for the past 2000 yr. Although Cox Island stratigraphy includes 9–12 peat-mud contacts much like those commonly inferred to record a series of megathrust earthquakes on the temperate coasts of this and other subduction zones, our mapping of stratigraphic contacts beneath tidal marshes near the river, lithologic descriptions of cores and outcrops, qualitative diatom analysis, quantitative foraminiferal analysis using a Bayesian transfer function, and 14C dating of the contacts failed to confirm the interpretation that any of the contacts of the past 2000 yr formed through sudden subsidence during great earthquakes.

Based on the youngest peat-mud contact’s (contact A) distinctness, wide (>400 m) distribution, shallow depth (∼0.5–0.6 m), and overlying probable tsunami-deposited sand bed, we correlated it with similar evidence for the 1700 CE earthquake and tsunami along much of the subduction zone. However, means of reconstructions of coseismic subsidence across the contact using the Bayesian foraminiferal transfer function range from 0.1 m to 0.5 m, suggesting that subsidence was below the threshold for distinguishing coseismic subsidence contacts from those of nonseismic origins (Shennan et al., 2016). The minimal changes in diatom assemblages across the contact are consistent with lower reconstruction values, as are geophysical models of the 1700 CE rupture (Wang et al., 2013; Gao et al., 2018; Wirth and Frankel, 2019).

However, for the other 11 peat-mud contacts mapped among our most thoroughly studied cores and river outcrop, their more limited stratigraphic extent and minimal changes in lithology, foraminifera, and(or) diatoms across them are insufficient to distinguish the contacts from those formed through small, gradual, or localized changes in tide levels resulting from changes in river-channel or estuary configuration during river floods, storm surges, and gradual sea-level rise. Although no data preclude any contacts from being synchronous with a megathrust earthquake, the evidence is equally consistent with the 11 contacts recording RSL changes below the 0.5 m detection threshold for distinguishing coseismic from nonseismic changes (e.g., Shennan et al., 2016). Many of our modeled age intervals for the 11 contacts overlap with similar age intervals having more distinct evidence of earthquake subsidence at sites to the north and south; for example, contact G’s interval overlaps with intervals for one of the greatest earthquakes of the past 2000 yr at many sites (Fig. 12). However, interval errors are too large and the time gaps between intervals too short to use interval overlap to infer an earthquake origin for any of the 11 Cox Island contacts. If our inference of tsunami deposition for beds overlying the fifth youngest contact (Db; 790–670 cal yr B.P.) is correct, it probably records a tsunami from a megathrust earthquake that produced little measurable subsidence at Cox Island.

However, in its failure to detect evidence sufficient to support an earthquake origin for any of the Cox Island peat-mud contacts, our study limits subsidence during megathrust earthquakes along this part of the subduction zone to <0.5 m for the past 2000 yr. Because all geophysical rupture models of slip heterogeneity along the Cascadia subduction zone rely on microfossil-based reconstructions of coseismic land-level change, a <0.5 m limitation could help limit the landward extent of modeled megathrust ruptures or suggest whether or not the locations of high-frequency subevent earthquakes remain stable over multiple earthquake cycles (Wang et al., 2013; Gao et al., 2018; Wirth and Frankel, 2019). Such a limitation is also consistent with offshore basin margin highs or the subduction of seamounts influencing the distribution and strength of slip patches on megathrust ruptures along this part of the subduction zone (Wells et al., 2003; Tréhu et al., 2012; Wirth and Frankel, 2019).

This work was supported by the Earthquake Hazards Program of the U.S. Geological Survey (USGS) and by U.S. National Science Foundation awards: 1419824 to Horton and 1419846 to Hawkes. Collection and analysis of diatom samples by Sawai were supported by the Geological Survey of Japan and the Japan Society for the Promotion of Science of postdoctoral fellowships for research abroad. Horton was also funded by the Singapore Ministry of Education Academic Research Fund MOE2018-T2–1–030, the National Research Foundation Singapore, and the Singapore Ministry of Education, under the Research Centres of Excellence initiative. The National Ocean Sciences Accelerator Mass Spectrometry facility (NOSAMS) at Woods Hole Oceanographic Institution (WHOI) supported the analysis of five 14C accelerator mass spectrometer (AMS) samples during Hawkes’ WHOI NOSAMS postdoctoral fellowship. We thank the Nature Conservancy for permission to study the stratigraphy beneath the marshes of Cox Island, and Eileen Hemphill-Haley for the format of Figures 9 and 10. Andrew Kemp (Tuffs University) provided key help in the field. Laura Brophy (2009) generously provided much unpublished data on plant communities and surveyed elevations on Cox Island. Jamie Delano (USGS, Golden) helped with graphics. This manuscript was much improved through reviews by Tina Dura (Virginia Polytechnic Institute and State University) and three anonymous reviewers. This work is Earth Observatory of Singapore contribution #332. This article is a contribution to PALSEA3 (Palaeo-Constraints on Sea-Level Rise) and International Geoscience Program (IGCP) Project 639, “Sea Level Change from Minutes to Millennia.” Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.

1Supplemental Material. Includes tables and imagery showing core and sampling locations; figures showing stratigraphy at additional sites and results of transfer function reconstructions of elevation using diatom floras from core S; tables of foraminiferal and diatom data; summaries of previous investigations; the tidal marsh setting of our study site; methods of measuring sampling elevations; explanation of variance added to radiocarbon age errors; and listing of code for OxCal radiocarbon age models. Please visit https://doi.org/10.1130/GEOS.S.13151009 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: David E. Fastovsky
Associate Editor: Jose Miguel Hurtado
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.