Incorporating temporal and spatial variability of salt-marsh foraminifera into sea-level reconstructions

Foraminifera from salt-marsh environments have been used extensively in quantitative relative sea-level re- constructions due to their strong relationship with tidal level. However, the influence of temporal and spatial variability of salt-marsh foraminifera on quantitative reconstructions remains unconstrained. Here, we con- ducted a monitoring study of foraminifera from four intertidal monitoring stations in New Jersey from high marsh environments over three years that included several extreme weather (temperature, precipitation, and storm surge) events. We sampled four replicates from each station seasonally (four times per year) for a total of 188 samples. The dead foraminiferal assemblages were separated into four site-specific assemblages. After ac- counting for systematic trends in changes in foraminifera over time among stations, the distribution of foraminiferal assemblages across monitoring stations explained ~87% of the remaining variation, while ~13% can be explained by temporal and/or spatial variability among the replicate samples. We applied a Bayesian transfer function to estimate the elevation of the four monitoring stations. All samples from each station predicted an elevation estimate within a 95% uncertainty interval consistent with the observed elevation of that station. Combining samples into replicate- and seasonal-aggregate datasets decreased elevation estimate uncertainty, with the greatest decrease in aggregate datasets from Fall and Winter. Information about the temporal and spatial variability of modern foraminiferal distributions was formally incorporated into the Bayesian transfer function through informative foraminifera variability priors and was applied to a Common Era relative sea-level record in New Jersey. The average difference in paleomarsh elevation estimates and uncertainties using an informative vs uninformative prior was minimal (< 0.01 m and 0.01 m, respectively). The dead foraminiferal assemblages remained consistent on temporal and small spatial scales, even during extreme weather events. Therefore, even when accounting for variability of modern foraminifera, foraminiferal-based relative sea-level reconstructions from high marsh environments remain robust and reproducible.

fossil assemblages from sediment cores, commonly from high marsh sedimentary sequences, to produce continuous records of RSL at decimeter vertical resolution (e.g. Horton et al., 1999;Gehrels, 2000;Horton and Edwards, 2006;Kemp and Telford, 2015;Cahill et al., 2016). A detailed understanding of modern salt-marsh foraminifera distributions (e.g. de Rijk, 1995) is a necessary prerequisite for the development of a transfer function to reconstruct RSL. Salt-marsh foraminifera are typically described in the modern environment only on one occasion in time without replicate sampling (e.g. Scott and Medioli, 1978;Horton, 1999;Kemp et al., 2017b). However, salt-marsh foraminifera have been found to vary temporally on seasonal and interannual timescales (Buzas et al., 2002;Hippensteel et al., 2002;Martin et al., 2002;Horton and Edwards, 2003;Horton and Murray, 2006;Berkeley et al., 2008), as well as spatially at small (sub-meter) scales (Buzas, 1968;Swallow, 2000;Morvan et al., 2006;Kemp et al., 2011). These variations have been shown to affect elevation boundaries of foraminiferal zones by as much as 15% of the tidal range (e.g. Horton and Edwards, 2003). However, such temporal and spatial variability has not been formally quantified in transfer functions used for RSL reconstructions.
We conducted a three-year monitoring study of foraminifera from four intertidal stations of differing salinity regimes within the high marsh to assess seasonal and interannual changes, as well as small-scale spatial variability through replicate samples at each station. This experiment is one of the longest seasonal/interannual monitoring studies of salt-marsh foraminifera (e.g. Hippensteel et al., 2002;Horton and Edwards, 2003;Horton and Murray, 2006) and importantly includes extreme weather events (temperature, precipitation, and storm surge). First, we analyzed the foraminiferal data to estimate: (a) the variation in abundance across monitoring stations for each species; and (b) the proportion of variation explained by the monitoring stations and replicate samples. Second, we applied a Bayesian transfer function on all of the data for each monitoring station to obtain a station elevation estimate and to examine variation in elevation estimates over time and across replicates at each station. Third, we utilized the information about the combined spatial and temporal variability by employing the species variance estimates from our analysis to inform prior distributions in the Bayesian transfer function. These variance estimates are incorporated into a RSL reconstruction from southern New Jersey (Kemp et al., 2013) to examine differences in the RSL reconstruction with informative vs uninformative foraminifera variability priors.

Regional setting
The field study sites are located in an intertidal environment of southern New Jersey on the U.S. mid-Atlantic coast (Fig. 1). The southern New Jersey coast is characterized by a barrier island and lagoon system adjacent to the Atlantic Ocean. Inlets between the barrier islands and lagoons allow water exchange between the ocean and bays (Ferland, 1990).
Modern salt marshes with tidal channels form extensive gently sloping (< 1:1000) platforms along the coast of southern New Jersey (Ferland, 1990). Low marsh environments are dominated by Spartina alterniflora (tall form), while high marsh environments are dominated by Spartina patens, Spartina alterniflora (short form), and Distichlis spicata (Daddario, 1961). The brackish environment between the high marsh and freshwater upland is vegetated by Phragmites australis and Iva fructescens (Daddario, 1961;Stuckey and Gould, 2000). Our high marsh field sites are located near the Rutgers University Marine Field Station (Tuckerton, New Jersey), in the Mullica River-Great Bay estuary, which is one of the most pristine estuaries on the U.S. mid-Atlantic coast with minimal human disturbance due to a lack of agricultural and industrial development and low population density (Kennish, 2004). The 1474 km 2 watershed is part of the Jacques Cousteau National Estuarine Research Reserve and drains the Pinelands National Reserve (Kennish, 2004).
The southern New Jersey coast has semidiurnal tides with microtidal (< 2 m) ranges that vary between the ocean and lagoon side of the barriers. The tidal range in the Mullica River-Great Bay estuary that influences our field sites varies from 0.7 m (in Little Egg Harbor) to 1.1 m (near the mouth of Great Bay). Water exchange primarily occurs between the Atlantic Ocean and Little Egg Inlet leading into Great Bay (Chant et al., 2000).
Meteorological data for the region during the sampling timeframe of this study was obtained from the Jacques Cousteau National Estuarine Research Reserve meteorological station at Nacote Creek (NOAA National Estuarine Research Reserve System (NERRS), 2019), ~12 km from the monitoring stations (Fig. 2). Average monthly air temperatures over the three-year study period ranged from −5 to 25 °C, with lows each year in January and February and highs in July and August. The study period contained six months that were the statewide warmest on record in New Jersey from 1895 to 2018 (Office of the New Jersey State Climatologist): May, November, and December in 2015; August in 2016; and February and April in 2017. Total monthly precipitation ranged from < 20 mm to > 200 mm. In 2015, the statewide third driest May and the fourth wettest June on record in New Jersey from 1895 to 2018 were observed (Office of the New Jersey State Climatologist). The Jacques Cousteau National Estuarine Research Reserve water quality station in Great Bay (~2 km from Station 1; ~10 km from Stations 2, 3, and 4) provided sea surface temperatures and salinity. Average monthly sea surface temperatures ranged from 5 to 25 °C, exhibiting comparable annual fluctuations with lows each year in January and highs in July and August, and average monthly salinity ranged from 28 to 32 ppt. Tide gauge data from the Atlantic City tide gauge (~19 km from Station 1; ~27 km from Stations 2, 3, and 4) were obtained from the Permanent Service for Mean Sea Level (Holgate et al., 2013). Monthly mean sea level heights exhibited annual lows in February and March and annual highs in September and October due to natural intra-annual variability driven by the annual warming/cooling cycle of changing seasons as well as fluctuations in salinities, winds, and currents. A significant winter storm flooding event was recorded at Atlantic City in January 2016, which was the fourth highest historic crest for the tide gauge, according to the NOAA National Weather Service.

Sampling design
We established four monitoring stations from high marsh/high marsh-upland transition sites above Mean Higher High Water (MHHW) along a salinity gradient in the Mullica River-Great Bay estuary (Fig. 1). We chose to investigate high marsh intertidal sites because sea-level studies use high marsh sedimentary sequences where foraminiferal zones are narrower compared to the low marsh, providing more precise elevation estimates (e.g. Gehrels, 2000;Kemp et al., 2011). We measured porewater salinity in the upper 10 cm of the marsh surface at the four monitoring station plots at each sampling period using a handheld YSI meter (Table 1). Station 1 had the highest salinity, while Stations 2 and 3 had intermediate salinities, and Station 4 had the lowest salinity.
At each station, we sampled a 1 m × 1 m plot four times per year (September, December, March, June) from September 2014 to June 2017 to examine temporal variability of salt-marsh foraminifera. Station 4 was established in March 2015. Samples for foraminiferal analysis were of a standardized volume of 10 cm 3 (10 cm 2 by 1 cm thick) to allow comparison with similar studies (e.g. Scott and Medioli, 1980;Horton and Edwards, 2006;Kemp et al., 2012). Four replicate surface sediment samples were collected so that small-scale spatial variability could be assessed. A different quadrant of each plot was sampled during each sampling period, following Horton et al. (2017), so that each quadrant was only sampled once per year to allow recovery of the marsh surface. In addition, we sampled for foraminifera at all four stations two weeks after the significant winter storm flooding event J.S. Walker, et al. Marine Geology 429 (2020) 106293 in January 2016. We surveyed each sampling station to NOAA tidal benchmarks using a total station, where elevations were referenced to the North American Vertical Datum (NAVD88). We took multiple elevation measurements within the 1 m × 1 m plot at Station 3 due to the uneven surface topography at this site. The elevation at Station 1 was converted to tidal datum levels using VDatum and the NOAA-operated tide gauge at Great Bay, Shooting Thorofare (station number 8534319) located < 500 m from the station. To convert elevations from Stations 2, 3 and 4 to tidal datum levels, we deployed two automatic water-level loggers (Solinst Levelogger Edge) in tidal channels within 100 m of the stations and leveled them to NOAA tidal benchmarks. We correlated the water-level logger data with those recorded by the NOAA-operated tide gauge at Tuckerton Creek (station number 8534080) located ~1 km north of Stations 2, 3, and 4. Due to differences in tidal range between Station 1 and Stations 2, 3, and 4, we converted the tidal elevations into a standardized water level index (SWLI), following the approach of Horton et al. (1999).

Foraminiferal analysis
We counted live and dead foraminifera from four replicate samples at each monitoring station from each sampling period for all three years. We stained the modern foraminifera samples with rose Bengal immediately after collection to distinguish live and dead foraminiferal tests (Walton, 1952). Although in some cases rose Bengal may stain dead tests (e.g. Walker et al., 1974;Bernhard, 1988), it remains a generally reliable method for identifying live tests and is unlikely to Creek and water quality station in Great Bay, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2.
Meteorological data for the region during the three-year sampling timeframe of this study from the Jacques Cousteau National Estuarine Research Reserve meteorological station at Nacote Creek and water quality station in Great Bay showing air temperature, total monthly precipitation, sea surface temperature, and salinity, including periods of extreme weather events. Tide gauge data were obtained from the Atlantic City tide gauge through the Permanent Service for Mean Sea Level (Holgate et al., 2013) and exhibit a cyclical annual pattern.  Walker, et al. Marine Geology 429 (2020) 106293 affect the interpretation of dead assemblages (Murray and Bowser, 2000). We stored samples in a buffered ethanol solution and refrigerated them (Scott et al., 2001). Samples were wet sieved to isolate the 63-500 μm size fraction and then split into eight equal aliquots using a wet splitter (Scott and Hermelin, 1993). We counted foraminifera under a binocular microscope while immersed in distilled water. Identifications of foraminifera were confirmed with type specimens at the National Museum of Natural History, Smithsonian Institute, Washington, D.C. Plate 1 shows scanning electron microscope (SEM) images of foraminifera from this study. We grouped specimens of the genera Haplophragmoides and of the genera Ammobaculites due to difficulties in identifying these genera to the species level (Kemp et al., 2009). Although live specimens were counted, only the dead assemblages are used in our analyses. Modern dead assemblages (compared to live or live plus dead assemblages) have been used for sea-level studies, because they most resemble subsurface assemblages and they are thought to minimize temporal variability in modern distributions (e.g. Horton, 1999;Horton and Edwards, 2003;Morvan et al., 2006).
Foraminiferal abundances are quoted as the mean with a 1σ standard deviation to quantify the variability of the data around the mean. All dead foraminifera data and Coefficients of Variation (CV; the ratio of the standard deviation to the mean) to measure the variability among replicate samples are presented in Appendix A.

Statistical analysis
We used partitioning around medoids (PAM) to analyze the composition of foraminifera assemblages present at the four monitoring stations (Kaufman and Rousseeuw, 1990). We applied PAM to the entire dataset of raw counts of dead foraminifera, using four groups to correspond to the four monitoring stations, and graphically represented the data with a silhouette plot (Rousseeuw, 1987). Analyses were completed using the 'cluster' package in R. Silhouette widths between −1 and 1 provide an estimate of a sample's classification. Values close to 1 indicate that the sample was assigned to an appropriate group where within group dissimilarity was less than the dissimilarity among J.S. Walker, et al. Marine Geology 429 (2020) 106293 the groups. Values close to −1 indicate that the sample was not appropriately classified. We incorporated variance estimates related to the spatial and temporal variation of foraminifera into a Bayesian transfer function. We first produced a regional dataset to use in the transfer function by merging our data with a modern foraminifera dataset from southern New Jersey (Kemp et al., 2013) and a new unpublished modern foraminifera dataset from Cheesequake State Park in northern New Jersey (Walker, 2019). The taxonomy was harmonized by combining Jadammina macrescens and Balticammina pseudomacrescens, and Trochammina inflata and Siphotrocha lobata. These species are often combined in sea-level transfer functions to avoid inconsistencies in taxonomic identifications (Kemp et al., 2018). Similar to PAM, the input data to the Bayesian transfer function was raw counts of dead foraminifera.
We performed three different analyses. First, we carried out a species variance analysis to identify the variability of each foraminiferal species across monitoring stations and among replicate samples. This was achieved by using the raw foraminifera counts (all samples over all three years) from all four of our study sites to estimate the overall species variance for each foraminiferal species. The species variance analysis was set up to: a) capture the underlying trends exhibited by each species across all monitoring stations over time; b) capture how these trends vary among each individual station; and c) capture the residual variation or noise outside of any underlying trend that may be present over time due to the replicate samples taken at each sampling period.
Second, we applied a Bayesian transfer function that employs foraminifera and a secondary proxy, bulk sediment stable carbon isotopes (Cahill et al., 2016;Kemp et al., 2017b). We used the regional dataset with the foraminifera data (raw counts) from each sampling period from each monitoring station, including replicates, to provide an elevation estimate for each station and to examine variability in these elevation estimates over the three years. Additionally, replicate sample assemblages were combined by adding together counts of all taxa to produce a replicate-aggregate elevation estimate for each sampling period for each station to analyze the effect of small-scale spatial variability. All sample assemblages for each season within each station were also combined across all three years to produce a seasonal-aggregate elevation estimate for each season for each station to analyze seasonal variability.
Here we provide an overview of the Bayesian transfer function. We outline notation for the data as follows: • y m are observed modern foraminifera abundances. y il m is the abundance of species l in surface sample i, with i = 1, …, N • T i m is the total of the species abundances for surface sample i. • e m are the observed standardized modern elevations. e i m is the elevation for surface sample i.
• y f are observed fossil foraminifera abundances. y jl f is the abundance of species l in fossil sample j, with j = 1, …, M • T j f is the total of the species abundances for fossil sample j.
A multinomial likelihood is assumed for the modern species abundances y il m as follows: where p il is the probability of finding species l at elevation i. The probability parameters p il m are estimated as a function of a latent response λ jl : where f is a softmax transformation used to preserve the sum and boundary constraints of the probabilities for each sample i. The latent response vector λ l contains the response for species l across all samples and is modeled as a function of elevation: where g l is a P-spline (e.g. de Boor, 1978;Dierckx, 1993) function that governs the shape of the response curve of species l. The error term ϵ l is added to the P-spline for each species to account for over/under dispersion in the raw data and σ l 2 is a species-specific variance term.
Third, we incorporated the overall species variance estimates and uncertainty from the species variance analysis into the Bayesian transfer function by providing informative priors for the species-specific variance terms (σ l 2 ) (subsequently referred to as the foraminifera variability prior) when reconstructing RSL. Specifically, a truncated tdistribution prior is placed on σ l such that, for the l th species The hyperparameters θ l and η l 2 control the magnitude and the uncertainty of σ l and the degrees of freedom ϑ = 1. The species variance analysis provided informed estimates for these hyperparameters.
To illustrate the influence of additional information regarding foraminiferal variability, we compared the transfer function results using informative foraminifera variability priors with the original Bayesian transfer functions of Cahill et al. (2016) and Kemp et al. (2017b) that had uninformative foraminifera variability priors.
We applied the Bayesian transfer function with informative and uninformative foraminifera variability priors to a RSL record from Kemp et al. (2013) from southern New Jersey (Fig. 1). The Bayesian transfer function produced SWLI estimates of paleomarsh elevation (PME) using fossil foraminifera abundances (raw counts) from a sediment core. The same modeling set up is assumed for fossil abundances, using the f subscript to refer to fossil data and parameters, as follows: The fossil elevations are contained within the vector e f , which has a prior distribution: where e j f is the fossil elevation for core sample j. The mean of the prior distribution for e j f has a uniform prior. In New Jersey, a secondary proxy is available (bulk sediment stable carbon isotopes) that can provide constraints on the elevational range of the fossil samples, which was used in the RSL record from Kemp et al. (2013). Therefore, a j and b j are fixed at the minimum and maximum elevations suggested by the secondary proxy for sample j. We converted SWLI estimates from the Bayesian transfer function to meters relative to Mean Tide Level (MTL). The PME estimates were subtracted from their sample altitude to obtain RSL. Finally, when combined with sample ages, we produced a probabilistic RSL reconstruction using an Errors-In-Variables Integrated Gaussian Process (EIV-IGP) model (Cahill et al., 2015) that accounts for the vertical and chronological uncertainties of the RSL data. We examined differences in past RSL and rates of RSL change in the southern New Jersey reconstruction with informative vs uninformative foraminifera variability priors.

Foraminiferal distributions
Over the three-year study, 72,804 modern foraminiferal tests were enumerated from the four stations, consisting of 14 agglutinated species from 188 samples, including the post-storm flooding samples. Across the stations, 50 samples were taken through time with 138 spatial replicate samples through the sampling period. Sample total count sizes ranged from 58 to 1196 tests/10 cm 3 with an average of 389 ± 221 tests/10 cm 3 (1σ) (Fig. 3). The total count of foraminifera was greatest in Year 3, when 26,194 tests were enumerated (average 409 ± 200 tests/10 cm 3 ). Furthermore, the maximum total counts of foraminifera were found in the last sampling date, Summer of Year 3, when a total of 7552 dead foraminifera were identified from the four replicates. The average range in count size among replicates was 235 ± 150 tests.
We used PAM analysis with four groups to examine the composition of the entire foraminifera dataset, including the post-storm flooding samples. The average silhouette width is 0.56, meaning the samples fit  Partitioning around medoids (PAM) analysis with four groups, showing four site-specific foraminiferal assemblages. All but 12 samples were assigned to a group corresponding to the monitoring station they were sampled from. J.S. Walker, et al. Marine Geology 429 (2020) 106293 well into four groups because the value is close to 1 (Fig. 4). The samples from Station 1 and Station 4 each fit into a group with an average silhouette width of 0.62, and Station 2 and Station 3 samples each fit into a group with an average silhouette width of 0.51. The higher silhouette widths for Stations 1 and 4 show that the samples in those groups are more similar to each other, and therefore their foraminifera assemblages are more consistent over time and space. All but 12 of 188 samples were assigned to a group corresponding to the monitoring station they were sampled from, indicating that each station has a site-specific assemblage. These 12 samples all had low silhouette widths (< 0.35) and did not as clearly belong to any one group. 10 of these 12 samples were from Stations 2 and 3, which had the groups with lower silhouette widths, and therefore these stations have less consistent assemblages over time and space. We note that the poststorm flooding samples were assigned to a group corresponding to the station they were sampled from.

Station 1 foraminiferal variability
Station 1 is a high marsh, high salinity site primarily vegetated by Spartina alterniflora (short form) with an average salinity over the threeyear sampling period of 39.6 ± 8.8 psu (1σ) and an elevation of 212-223 SWLI units (0.68 ± 0.03 m MTL, 1σ). Station 1's total count size ranged from 58 to 616 tests/10 cm 3 with an average of 229 ± 112 tests/10 cm 3 (1σ) (Fig. 5). The sample after flooding in January 2016 had a total count of 110 tests/10 cm 3 . The total count of foraminifera was greatest in Year 3 when 4618 tests were enumerated (average 289 ± 148 tests/10 cm 3 ). Furthermore, the maximum abundance of foraminifera was found in the last sampling date, Summer of Year 3, when 1596 dead foraminifera were identified from the four replicates.
Year 2 had the lowest (96 ± 37 tests) and Year 3 had the highest (284 ± 180 tests) average range in count size among replicates. The maximum range of foraminifera among replicates was found in Spring of Year 3 with 68 to 528 tests/10 cm 3 .
We identified 11 foraminiferal species from Station 1 that were dominated by T. inflata, J. macrescens, and T. comprimata, including for the sample after flooding in January 2016 (Fig. 5). T. inflata was the dominant species in 47 of the 48 samples from Station 1. T. inflata's total count varied from 22 to 436 tests/10 cm 3 (average 144 ± 80 tests/10 cm 3 ). J. macrescens had a total count that varied from 6 to 166 tests/10 cm 3 (average 53 ± 34 tests/10 cm 3 ). Similar to all of the species combined, both T. inflata and J. macrescens' total counts were greatest in Year 3 (3080 and 1148 total tests, respectively), and neither species exhibited a seasonal pattern. T. comprimata had a total count that varied from 2 to 132 tests/10 cm 3 (average 23 ± 25 tests/ 10 cm 3 ). The annual T. comprimata total count remained relatively stable from 400 total tests in Year 1 (average 25 ± 29 tests/10 cm 3 ) to 362 total tests in Year 3 (average 23 ± 31 tests/10 cm 3 ). T. comprimata did, however, exhibit a seasonal pattern. The maximum total counts of each year were found in Spring when 226 (2015), 124 (2016), and 192 (2017) tests were identified from the four replicates. The replicate sample ranges for T. comprimata were also greatest in the Spring sampling periods for all three years.

Station 2 foraminiferal variability
Station 2 is primarily vegetated by Spartina patens and had an average salinity over the three-year sampling period of 13.7 ± 4.9 psu and an elevation of 211-224 SWLI units (0.68 ± 0.03 m MTL). Station 2's count size ranged from 186 to 1130 tests/10 cm 3 with an average of 500 ± 225 tests/10 cm 3 (Fig. 6). The sample after flooding in January 2016 had a total count of 338 tests/10 cm 3 . The total count of foraminifera was greatest in Year 1 (average 601 ± 269 tests/10 cm 3 ) with the maximum abundance of foraminifera in the first sampling Year 1 had the highest average range in count size (393 ± 190 tests). The greatest range of foraminifera among replicates was found in the first sampling date, Fall of Year 1, with 478 to 1130 tests/10 cm 3 .
We identified 10 foraminiferal species from Station 2 that were dominated by J. macrescens, B. pseudomacrescens, and T. comprimata, including for the sample after flooding in January 2016 (Fig. 6). J. macrescens was the dominant species in 35 of the 48 samples. J. macrescens' total count varied from 46 to 656 tests/10 cm 3 (average 239 ± 141 tests/10 cm 3 ) and did not exhibit a seasonal pattern. B. pseudomacrescens had a total count that varied from 14 to 456 tests/ 10 cm 3 (average 157 ± 91 tests/10 cm 3 ). B. pseudomacrescens' total count was relatively stable seasonally and annually except in Spring of Year 1 when 1296 tests were identified from the four replicates. T. comprimata had a total count that varied from 8 to 346 tests/10 cm 3 (average 76 ± 64 tests/10 cm 3 ) and did not exhibit a seasonal pattern. Similar to all of the species combined, T. comprimata's total count was greatest in Year 1 with the maximum abundance found in the first sampling date, Fall of Year 1.

Station 3 foraminiferal variability
Station 3 is primarily vegetated by Spartina patens and Distichlis spicata, but borders a Phragmites australis flora, and had an average salinity slightly lower than Station 2 of 13.1 ± 4.9 psu. Station 3 has an uneven surface topography and its elevation ranges from 192 to 227 SWLI units (0.46-0.60 m MTL). Station 3 had count sizes ranging from 212 to 1196 tests/10 cm 3 with an average of 553 ± 198 tests/10 cm 3 (Fig. 7). The sample after flooding in January 2016 had a total count of 354 tests/10 cm 3 . The total count of foraminifera was greatest in Year 2 (average 622 ± 178 tests/10 cm 3 ) with the maximum abundance of foraminifera in Fall of Year 2 (3188 dead foraminifera from the four replicates). Of the four stations, Station 3 had the highest average range in count size among replicates of 339 ± 157 tests. Year 1 had the highest average range in count size among replicates at 407 ± 235 tests with the greatest range of foraminifera among replicates in the first sampling date, Fall of Year 1 (504 to 1196 tests/10 cm 3 ).
We identified 12 foraminifera species from Station 3 that were dominated by T. comprimata, B. pseudomacrescens, and Haplophragmoides spp. (Fig. 7). The post-flooding sample was dominated by T. comprimata, Haplophragmoides spp., and T. inflata. T. inflata was the overall fourth dominant species over the sampling period at Station 3 and B. pseudomacrescens was present in the post-flooding sample. T. comprimata was the dominant species in 34 of the 48 samples. T. comprimata's total count varied from 60 to 380 tests/10 cm 3 (average 190 ± 74 tests/10 cm 3 ). Similar to Station 1, T. comprimata exhibited a seasonal pattern; however, an opposite pattern was observed. Maximum abundances each year were found in Fall and Winter, and minimum abundances each year were found in Spring and Summer. B. pseudomacrescens had a total count that varied from 38 to 302 tests/ 10 cm 3 (average 128 ± 60 tests/10 cm 3 ). B. pseudomacrescens' total count was relatively stable interannually except in Fall of Year 2 when a total of 910 tests were identified from the four replicates. Seasonally, the total count was highest in the Fall for all three years. Haplophragmoides spp. had a total count that varied from 12 to 266 tests/ 10 cm 3 (average 92 ± 60 tests/10 cm 3 ) and did not exhibit a seasonal

Station 4 foraminiferal variability
Station 4 is a high marsh-upland transition, brackish site vegetated by Phragmites australis and had the lowest salinity with an average over the three-year sampling period of 2.9 ± 2.0 psu and the highest elevation of the four monitoring stations at 248-264 SWLI units (0.61 ± 0.03 m MTL). Station 4 was established in March 2015; therefore, Year 1 refers to data only from Spring and Summer 2015. Station 4 had count sizes ranging from 88 to 456 tests/10 cm 3 with an average of 252 ± 86 tests/10 cm 3 (Fig. 8). The sample after flooding in January 2016 had a total count of 344 tests/10 cm 3 . The total count of foraminifera remained stable through all three years from an average of 265 ± 100 tests/10 cm 3 in Year 1 to an average of 242 ± 104 tests/ 10 cm 3 in Year 2 to an average of 256 ± 60 tests/10 cm 3 in Year 3. Of the four monitoring stations, Station 4 had the smallest average range in count size among replicates of 156 ± 68 tests. Year 1 had the highest average range in count size among replicates (245 ± 7 tests).
We identified 10 foraminifera species from Station 4 that were dominated by J. macrescens, A. inepta, B. pseudomacrescens, and Haplophragmoides spp. (Fig. 8). The dominant species in the postflooding sample were J. macrescens, A. inepta, Haplophragmoides spp., and M. petila. M. petila was the overall fifth dominant species over the sampling period at Station 4 and B. pseudomacrescens was present in the post-flooding sample. J. macrescens was the dominant species in all 40 samples. J. macrescens' total count varied from 58 to 294 tests/10 cm 3 (average 166 ± 64 tests/10 cm 3 ) and did not exhibit a seasonal pattern. The total count was lowest in the Fall for both years that included samples in September. A. inepta had a total count that varied from 2 to 124 tests/10 cm 3 (average 30 ± 30 tests/10 cm 3 ) and did not exhibit a seasonal pattern. B. pseudomacrescens and Haplophragmoides spp. had total counts that varied from 0 to 56 tests/10 cm 3 (average of 15 ± 14 tests/10 cm 3 ). Haplophragmoides spp.'s total count showed annual increases from an average of 7 ± 5 tests/10 cm 3 in Year 1 to 22 ± 10 tests/10 cm 3 in Year 3. In contrast, B. pseudomacrescens' total count showed annual decreases during the study period from an average of 24 ± 14 tests/10 cm 3 in Year 1 to 6 ± 5 tests/10 cm 3 in Year 3. Neither B. pseudomacrescens nor Haplophragmoides spp. exhibited a seasonal pattern.

Bayesian transfer function elevation estimates
We used a Bayesian transfer function to estimate the elevation of each monitoring station from our foraminiferal data for each sampling period at each monitoring station. We illustrate the analysis for Station 1 in Fig. 9, with the remaining stations summarized in Table 2 and Supplementary Fig. 1, 2, and 3. Under a 95% uncertainty interval, all samples from each monitoring station, including the samples after flooding in January 2016, predicted a SWLI estimate within the observed elevation range of that station. Station 1 had the smallest range (197-219) of SWLI estimates and smallest average SWLI uncertainty (26.5) among the stations, and Station 3 had the largest range (133-240) of SWLI estimates and largest average SWLI uncertainty (36.5) among the stations.
SWLI estimates have a relationship with count size and the presence or greater number of rare species. Samples with a lower count size compared to replicates taken at the same time often have an anomalous SWLI estimate. For example, in Summer of Year 1 at Station 2 (observed elevation = 211-224), the four replicate samples had comparable foraminifera assemblages, but count sizes of 480, 554, 228, and 488 tests and SWLI estimates of 226, 222, 195, and 222, respectively. The presence of rare species, even in the smallest numbers, such as M. fusca or Ammobaculites spp. decreases a sample's SWLI estimate, while M. petila increases a sample's SWLI estimate. For example, in Spring of Year 2 at Station 1 (observed elevation = 212-223), the four replicate samples had SWLI estimates of 212, 211, 213, and 201. The replicates had comparable assemblages, except the fourth sample had the presence of M. fusca (2 tests).
Replicate sample assemblages were added together to produce a replicate-aggregate SWLI estimate for each sampling period for each station to examine the influence of combining replicate samples taken from a small-scale spatial area on elevation estimates. The range in all SWLI estimates decreased in all four stations when using the replicateaggregate dataset with the largest decrease at Station 4 (74 SWLI units). The replicate-aggregate dataset decreased the average uncertainty in the SWLI estimates compared to using the full dataset at Stations 1, 2, and 4, with the greatest decrease of 11 SWLI units at Station 4. Of the four seasons, the average uncertainty in the SWLI estimates decreased when using the replicate-aggregate dataset compared to the full dataset for all seasons at Stations 1, 2, and 4, with the greatest decrease in Fall for Stations 1 and 2 and in Winter for Station 4. The replicate-aggregate uncertainty at Station 3 decreased only in Fall (by 1.5 SWLI units) and Summer (by 1 SWLI unit).
All sample assemblages for each season within each station were also added together to produce a seasonal-aggregate SWLI estimate for each season for each station to examine the influence of seasonal variability of foraminifera on elevation estimates. The average of the Fall, Winter, Spring, and Summer seasonal-aggregate SWLI estimates for each monitoring station did not significantly change from the average SWLI estimate using the full dataset or the replicate-aggregate dataset. The range in all SWLI estimates further decreased from the replicate-aggregate dataset at Stations 1, 3, and 4 when using the seasonal-aggregate dataset with the largest decrease at Station 3 (43 SWLI units). The Fall and Winter seasonal estimates for each monitoring station were within 3 SWLI units of the observed station elevations, while the Spring and Summer seasonal estimates were up to 15 SWLI units different from the observed station elevations. The seasonal-aggregate dataset also lowered the uncertainty in the SWLI estimates compared to using the full dataset or the replicate-aggregate dataset for all seasons except Summer at Stations 1 and 4. The uncertainty was reduced only in Winter (by 3.5 SWLI units) and Summer (by 1 SWLI unit) at Station 2 and only in Winter at Station 3 (by 4.5 SWLI units).

Informing variability in the Bayesian transfer function
The species variance analysis of the entire raw foraminifera dataset illustrates that the variation across monitoring stations made up 87% (95% credible interval of 56-95%) of the total variation in the foraminiferal dataset, while the remaining 13% (5-44%) of the variation can be explained by temporal and/or spatial variability among the replicates. The combination of the variation across monitoring stations and the variation among replicates contributed to an overall variation term that was estimated for the dominant species and subsequently incorporated into the Bayesian transfer function.
We incorporated the species-specific temporal and spatial uncertainty from the species variance analysis into the Bayesian transfer function by providing informative priors for the relevant variation parameters in the model (foraminifera variability prior). This is in addition to a bulk sediment stable carbon isotope prior (Cahill et al., 2016) used to inform elevation estimates. The Bayesian transfer function including informative/uninformative foraminifera variability priors was applied to the fossil foraminiferal data from a southern New Jersey RSL record of Kemp et al. (2013) (Fig. 10). The PME estimates from the transfer function with both the uninformative and informative foraminifera variability priors were consistent with one another. The average difference in PME estimates was < 0.01 m and all PME estimates overlapped within the 95% uncertainty interval. Furthermore, the average difference in PME estimate uncertainties was 0.01 m. The EIV-IGP model found very similar RSL change over the past ~1000 years: 1.66 m (95% credible interval of 1.40-1.89 m) rise with the uninformative foraminifera variability prior and a 1.62 m (95% credible interval of 1.36-1.86 m) rise with the informative prior (Fig. 10). Furthermore, the average difference in rate predictions was 0.04 ± 0.12 mm/yr and all rate predictions overlapped within the 95% uncertainty interval.

Foraminiferal distributions
The dead foraminiferal distributions from the high marsh and high marsh-upland transition monitoring stations in the Mullica River-Great Bay estuary are similar to other studies in New Jersey and on the U.S. Atlantic coast (e.g. Culver et al., 1996;Hippensteel et al., 2000;Kemp et al., 2009Kemp et al., , 2011. The four monitoring stations were dominated by T. inflata, J. macrescens, T. comprimata, B. pseudomacrescens, Haplophragmoides spp., and A. inepta. Although Stations 1, 2, and 3 are all located in high marsh above MHHW and Station 4 in high marsh-upland transition above MHHW, the foraminiferal assemblages and dominant foraminifera differ. The species variance estimates illustrated that the variation across monitoring stations made up ~87% of the total variation in the foraminiferal dataset, demonstrating unique site-specific assemblages. High marsh assemblages of dead foraminifera have been shown to vary both among and within regions (e.g. Ellison and Nichols, 1976;Wright et al., 2011;Kemp et al., 2013). While salt-marsh foraminifera distributions are strongly linked with tidal elevation (e.g. Horton and Edwards, 2006;Kemp et al., 2013), variability in foraminiferal assemblages among high marsh sites may also be influenced by secondary environmental factors such as salinity (e.g. de Rijk, 1995;Nikitina et al., 2003;Kemp et al., 2009;Wright et al., 2011;Kemp et al., 2012Kemp et al., , 2013. Our four monitoring stations exhibit a salinity gradient, which Replicate-aggregate SWLI estimates are shown as red data points on top of the estimated SWLI using the full dataset. (C) Seasonal SWLI estimates are shown for each season (Fall, Winter, Spring, Summer) using the full dataset, a replicate-aggregate dataset, and a seasonal-aggregate dataset. (D) Similarly, seasonal SWLI uncertainties in elevation estimates are shown for each season using the three different datasets. Equivalent analysis for Stations 2, 3, and 4 can be found in Supplementary Fig. 1, 2, and 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) varies by up to 35 psu among the stations. Along the Atlantic coast of North America, Wright et al. (2011) demonstrated the spatial differences in high marsh foraminifera assemblages among and within regions. For example, in the Chesapeake Bay region, middle and highmarsh environments consist of assemblages of variable proportions of several dominant species, which were correlated with salinity gradients (Ellison et al., 1965;Ellison and Nichols, 1976). In North Carolina and New Jersey, Kemp et al. (2009Kemp et al. ( , 2013 also found sub-regional groups of foraminifera from high marsh environments where spatial differences in species composition likely reflected the distribution of salinity regimes of the region due to the balance between marine tidal influence and freshwater input from rivers at individual sites. Additionally, variability among marsh assemblages could be linked to morphotypes of individual species that correspond to factors such as salinity (Scott and Medioli, 1980). Station 1, located in Spartina alterniflora (short form) high marsh, has the highest salinity (39.6 ± 8.8 psu) of the four monitoring stations and a foraminiferal assemblage dominated by agglutinated species T. inflata, J. macrescens, and T. comprimata. Similar foraminiferal assemblages have been observed in high-salinity, high-marsh environments. On the U.S. Atlantic coast, T. inflata has been recognized as a dominant species in the middle and high marsh (e.g. Hippensteel et al., 2000;Kemp et al., 2009;Kemp et al., 2012). In North Carolina, Kemp et al. (2009) also found T. inflata-dominated assemblages associated with higher salinity sites. Kemp et al. (2012) recognized distinct high marsh groups with varying proportions of T. inflata, J. macrescens, and T. comprimata in New Jersey.
Stations 2 and 3 are both located in high marsh environments with comparable, moderate salinities (Station 2 13.7 ± 4.9 psu, and Station 3 13.1 ± 4.9 psu), but have differing assemblages. Station 2 has a foraminiferal assemblage dominated by agglutinated J. macrescens, B. pseudomacrescens, and T. comprimata, and Station 3 has a foraminiferal assemblage dominated by agglutinated B. pseudomacrescens, T. comprimata, and Haplophragmoides spp. The marsh surface of Station 3 has an uneven surface topography, which results in a larger elevational range , compared to the elevation at Station 2 (211-224 SWLI units). De Rijk and Troelstra (1997) noted variability in assemblages due to microtopography with pond holes changing salinities. Further, 10 of the 12 samples that were not assigned to a PAM analysis group corresponding to the monitoring station they were sampled from were from Stations 2 and 3. This incorrect grouping and the lower silhouette widths of these 10 samples suggests a less consistent assemblage over time.
Similar foraminiferal assemblages have been observed in high marsh/low salinity salt marshes that correspond to both Stations 2 and 3. In New Jersey and North Carolina, high marsh assemblages have been dominated by J. macrescens, T. comprimata, and Haplophragmoides spp. (Kemp et al., 2009(Kemp et al., , 2012(Kemp et al., , 2013. Assemblages dominated by T. comprimata, which is a dominant species at both Stations 2 and 3, have been associated with lower salinity sites (Kemp et al., 2013). Kemp et al. (2012) found lower abundances of B. pseudomacrescens in the high marsh in New Jersey, although it is a dominant species at both Stations 2 and 3. B. pseudomacrescens has not been recorded in North Carolina or Virginia (Spencer, 2000;Kemp et al., 2009;Wright et al., 2011), but has been more prevalent in New England and Newfoundland, Canada (de Rijk, 1995;de Rijk and Troelstra, 1997;Gehrels and van de Plassche, 1999;Edwards et al., 2004;Wright et al., 2011). In New Jersey, Haplophragmoides spp., which is a dominant species at Station 3, has been found to be a dominant species in high marsh and transitional high marsh-upland environments, often above MHHW (Kemp et al., 2012(Kemp et al., , 2013, and in Massachusetts, Haplophragmoides spp. has been associated with low salinity, high marsh settings (de Rijk, 1995;de Rijk and Troelstra, 1997).
Station 4, located in a Phragmites australis high marsh-upland transition zone, has the lowest salinity (2.9 ± 2.0 psu) and highest elevation (248-264 SWLI) of the four monitoring stations. The foraminiferal assemblage is dominated by agglutinated species J. macrescens, which is consistent with other studies on the U.S. mid-Atlantic coast, which found maximum abundances of J. macrescens in high marsh-upland transition environments with low salinities (e.g. Spencer, 2000;Nikitina et al., 2003;Robinson and McBride, 2006;Horton and Culver, 2008;Kemp et al., 2009). Beginning with the work of Medioli (1978, 1980), salt-marsh foraminifera assemblages dominated by J. macrescens have been considered the highest elevational zone at the high marsh-upland transition. A. inepta, B. pseudomacrescens, and Haplophragmoides spp. are also found in high abundances at Station 4. In low salinity, brackish high marsh-upland transition environments, A. inepta has been found to be a dominant species in New Jersey and North Carolina (Scott et al., 2001;Culver and Horton, 2005;Kemp et al., 2009;Kemp et al., 2013). In New Jersey, Kemp et al. (2012Kemp et al. ( , 2013 found low abundances of B. pseudomacrescens in the highest marsh zones and found greater abundances of Haplophragmoides spp. in transitional environments above MHHW.

Temporal and spatial variability
Most studies of temporal and small-scale spatial variability of foraminiferal assemblages in salt marshes have focused on living populations (Lynts, 1966;Buzas, 1968Buzas, , 1970Schafer, 1971;Swallow, 2000;Buzas et al., 2002;Berkeley et al., 2008); few have considered dead assemblages Horton and Edwards, 2003;Horton and Murray, 2006;Morvan et al., 2006;Kemp et al., 2011). However, many studies have stressed the importance of replicate  Walker, et al. Marine Geology 429 (2020) 106293 sampling when studying modern foraminifera (Buzas, 1969;Schafer, 1971;Murray and Alve, 2000;Buzas et al., 2002). The dead foraminiferal assemblages and dominant species from the Mullica River-Great Bay estuary remained relatively consistent temporally at each monitoring station during the three-year study period. The variations in annual total count that were observed did not have a relationship with time. Station 1's total count was greatest in Year 3 whereas for Station 2 and 3 it was in Years 1 and 2, respectively. Seasonal variability of individual foraminiferal species was also inconsistent among the four monitoring stations. For example, although T. comprimata exhibited a seasonal pattern at two of the three monitoring stations where it was a dominant species, the pattern was not the same between stations. Live assemblages of foraminifera commonly show seasonal variability (e.g. Buzas and Hayek, 2000;Swallow, 2000). However, there was an absence of a translation from live foraminifera seasonal life cycles into the dead assemblage at Mullica River-Great Bay estuary, which has been documented elsewhere (Horton andMurray, 2006, 2007;Morvan et al., 2006). Indeed, other studies have also found differing densities of dead foraminiferal assemblages between years and no clear annual pattern in total counts (e.g. Hippensteel et al., 2002;Horton and Edwards, 2003).
Foraminiferal assemblages and dominant species from the Mullica River-Great Bay estuary also remained consistent on small spatial scales at each monitoring station over the study period. Kemp et al. (2011) also found that dead foraminifera in high marsh environments exhibited a non-patchy distribution. Total counts, however, varied among replicate samples. For example, the range of total count among replicates at Station 3 was 407 ± 235 tests/10 cm 3 . Observed smallscale spatial variability in total count may be due to a variety of environmental factors affecting the live assemblage, including response to predation (Buzas, 1978(Buzas, , 1982, reproduction (Stouff et al., 1999), availability of food resources (Alve and Murray, 2001;Fontanier et al., 2003) and species interactions (Buzas, 1968;Hayward et al., 1996;Scott et al., 2001), which then may be influencing the distribution of the dead assemblages. The length of this study also provided the opportunity to examine the influence of extreme weather events on foraminifera assemblages. Storm flooding events did not affect dead foraminiferal assemblages and dominant species. The marsh assemblages experienced no change or rapidly recovered two weeks after the winter flooding event in January 2016. There was no evidence for the influence of flooding on the marshes, such as the presence of overwash material, which can be deposited in marsh environments and identified using changes in foraminifera assemblages (e.g. Hippensteel and Martin, 1999; Culver et al., Fig. 10. Comparison of southern New Jersey relative sea-level Bayesian transfer function reconstruction (Kemp et al., 2013) with informative vs uninformative foraminiferal variability priors to account for temporal and spatial uncertainties in modern foraminiferal distributions. (A) Paleomarsh elevation (PME) estimates and uncertainties from the Bayesian transfer function are compared by core depth as a difference of estimates/uncertainties with the uninformative foraminiferal variability prior minus estimates/uncertainties with the informative prior. (B) An Errors-In-Variables Integrated Gaussian Process model (Cahill et al., 2015) compares the RSL record and rates of change through time with the uninformative vs informative foraminiferal variability prior.
2006; Pilarczyk et al., 2014). Dead foraminiferal counts did not show a strong correlation with extreme local climate events (Fig. 2), such as the record monthly air temperatures (May, November, December 2015; August 2016; February, April 2017) and precipitation extremes (Supplementary Fig. 4). For example, May 2015 was the statewide third driest May on record in New Jersey from 1895 to 2018, which was followed by the fourth wettest June on record the following month, and samples taken at the end of June showed no change in the dead assemblages. Further dead foraminifera did not show a strong correlation with intra-or interannual changes in sea surface salinity/temperature and mean sea level height (Supplementary Fig. 4).

Implications for sea-level studies
The temporal and spatial consistency, even after extreme weather events, of the dead foraminiferal assemblages in the Mullica River-Great Bay estuary is also reflected by the SWLI estimates for Stations 1-4 from the Bayesian transfer function ( Fig. 9; Supplementary Fig. 1,  2, and 3). The 95% uncertainty interval for each SWLI estimate contained the observed elevation for each sample from its corresponding monitoring station. However, the variability in total count and the presence of rare species of the foraminiferal assemblages influenced the elevation estimates and uncertainty from the Bayesian transfer function. Samples with a smaller count size compared to replicates taken at the same time often have a SWLI estimate that is anomalous compared to the other replicates, suggesting the importance of count size in quantitative studies of foraminifera. A simulation of the influence of count size showed a reduction of SWLI estimate uncertainties with increasing count sizes, which stabilizes with count sizes > 100 (Supplementary Fig. 5). Other studies also suggest counts of at least 100 are needed to fully capture non-dominant species within an assemblage (e.g. Buzas, 1990;Hayek and Buzas, 2010;Fatela and Taborda, 2002).
The presence or greater number of rare species (< 3 tests), especially of M. fusca or Ammobaculites spp., appears to consistently decrease a sample's SWLI estimate, while the presence or greater number of M. petila consistently increases a sample's SWLI estimate. These findings are consistent with the observations that M. fusca and Ammobaculites spp. are associated with lower elevations (e.g. Hippensteel et al., 2000;Edwards et al., 2004;Horton and Culver, 2008) and M. petila is associated with higher elevations (e.g. Scott and Medioli, 1978;Spencer, 2000;Kemp et al., 2009Kemp et al., , 2013. Therefore, rare species should be appropriately accounted for through a sufficient sample count size or through replicate sampling, especially when using raw counts rather than relative abundance data when all species are included in the analysis. Combining foraminiferal data from replicate samples decreased elevation estimate uncertainty, suggesting the addition of replicate samples provides a greater understanding of a modern site's foraminiferal distributions (e.g. Schafer, 1971;Murray and Alve, 2000;Buzas et al., 2002). Combining foraminiferal data from samples taken in the same seasons generally further lowered elevation estimate uncertainties compared to the replicate-aggregate dataset. Both replicateaggregate and seasonal-aggregate datasets had more accurate SWLI estimates, as well as lower uncertainties, in Fall and/or Winter compared to Spring or Summer. In a seasonal study of foraminifera, Horton and Edwards (2003) also found that the greatest transfer function precision was obtained using samples collected in the Winter months and the weakest in the Summer, and suggested a modern foraminifera dataset that includes samples spanning all seasons will provide the bestquality data for sea-level studies, which this study also supports.
For the first time, an informative prior has been developed to account for temporal and spatial variability of modern foraminifera to include in transfer functions to reconstruct RSL change. Incorporating a more informative foraminiferal variability prior into the Bayesian transfer function for the Kemp et al. (2013) southern New Jersey record, which uses a secondary proxy (bulk sediment stable carbon isotopes), minimally affected PME estimates due to the consistency of the Mullica River-Great Bay estuary modern foraminifera assemblages through time and space. All PME estimates, total RSL change and predicted rates of change with informative vs uninformative foraminiferal variability priors overlapped within a 95% uncertainty interval, indicating the minimal influence that temporal/spatial foraminiferal variability in high marsh environments has on RSL reconstructions. Therefore, although variability in foraminiferal assemblages has been shown to affect elevation boundaries of foraminiferal zones by as much as 15% of the tidal range (Horton and Edwards, 2003) and here we document spatial and temporal variations at high marsh locations, accounting for modern foraminiferal variability still provides consistent high marsh RSL reconstructions. High marsh sedimentary environments have been used for sea-level studies partly because of the consistency of high marsh foraminifera assemblages within individual sites and their narrow elevation zones compared to middle and low marsh environments (e.g. Gehrels, 2000;Kemp et al., 2011), which we demonstrate here through this monitoring study.

Conclusions
A detailed understanding of modern salt-marsh foraminiferal distributions is necessary to produce RSL reconstructions using foraminiferal-based transfer functions. We sampled four high marsh monitoring stations in the Mullica River-Great Bay estuary every three months over three years to examine seasonal and interannual changes and small-scale spatial variability in dead modern foraminiferal assemblages.
We found four site-specific assemblages where the variation across monitoring stations explained ~87% of the total variation in the foraminiferal dataset, while the remaining ~13% of the variation can be explained by temporal and/or spatial variability among replicate samples. Overall foraminiferal assemblages and dominant foraminiferal species at each monitoring station over the study period remained consistent both temporally and spatially among replicate samples, including after extreme weather events.
Combining replicate samples into a replicate-aggregate dataset lowered the uncertainty of elevation estimates and a seasonal-aggregate dataset further lowered elevation estimate uncertainties. Samples in the aggregate datasets from Fall and Winter months had more accurate elevation estimates, as well as lower uncertainties, compared to Spring and Summer.
Using a Bayesian transfer function with a modern foraminifera dataset for New Jersey, we found that under a 95% uncertainty interval, all samples from each monitoring station predicted a SWLI estimate within the observed elevation range of that station. Incorporating an informative foraminiferal variability prior to account for temporal and spatial changes in modern foraminiferal distributions into a RSL record in New Jersey, which also included a secondary proxy (bulk sediment stable carbon isotopes; Kemp et al., 2013), resulted in minimal changes in PME estimates and reconstructed RSL and rates of change.
RSL reconstructions using salt-marsh foraminifera rely on the consistency of foraminiferal assemblages in time and space in the modern environment. This study demonstrates that foraminiferal-based RSL reconstructions remain robust and reproducible even when accounting for temporal and spatial variability of salt-marsh foraminifera in the modern environment, including after extreme weather events. The informative foraminiferal variability prior could be applied to locations with similar high marsh foraminiferal assemblages to our study sites in southern New Jersey, such as elsewhere along the U.S. Atlantic coast, where modern foraminiferal assemblages may exhibit temporal and/or small-scale spatial variability.

Data availability
Datasets related to this article can be found online at https://rucore. libraries.rutgers.edu.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.