Drivers of 20th century sea‐level change in southern New Zealand determined from proxy and instrumental records

In this paper we present new proxy‐based sea‐level reconstructions for southern New Zealand spanning the last millennium. These palaeo sea‐level records usefully complement sparse Southern Hemisphere proxy and tide‐gauge sea‐level datasets and, in combination with instrumental observations, can test hypotheses about the drivers of 20th century global sea‐level change, including land‐based ice melt and regional sterodynamics. We develop sea‐level transfer functions from regional datasets of salt‐marsh foraminifera to establish a new proxy‐based sea‐level record at Mokomoko Inlet, at the southern tip of the South Island, and to improve the previously published sea‐level reconstruction at Pounawea, located about 110 km to the east. Chronologies are based on radiocarbon, radiocaesium, stable lead isotope and pollen analyses. Both records are in good agreement and show a rapid sea‐level rise in the first half of the 20th century that peaked in the 1940s. Previously reported discrepancies between proxy‐based sea‐level records and tide‐gauge records are partially reconciled by accounting for barystatic and sterodynamic components of regional sea‐level rise. We conclude that the rapid sea‐level rise during the mid‐20th century along the coast of southern New Zealand was primarily driven by regional thermal expansion and ocean dynamics.


Introduction
Late Holocene microfossil-based relative sea-level (RSL) reconstructions complement observational data from tide gauges and are crucial for understanding the rates and processes driving sea-level change from the pre-instrumental period to the present day (e.g. Gehrels et al., 2005Gehrels et al., , 2020Kemp et al., 2011;Barlow et al., 2013;Gehrels and Woodworth, 2013;Kopp et al., 2016;Walker et al., 2021). Southern Hemisphere proxy records (e.g. Frederikse et al., 2021) are rare in comparison to more numerous reconstructions from Europe and North America (e.g. Gehrels, 2000;Donnelly et al., 2004;Gehrels et al., 2005;Kemp et al., 2011Kemp et al., , 2018Long et al., 2014), which is perhaps surprising given that the Southern Hemisphere hosts most of the world's oceans (Dangendorf et al., 2019). Spatiotemporal patterns in these far-field sea-level data are also of critical importance for determining the processes driving Southern Hemisphere RSL change, including quantifying glacial massbalance fluctuations in response to Common Era climate change (Mitrovica et al., 2001;Gehrels and Woodworth, 2013;Gehrels et al., 2020).
Existing proxy reconstructions from Tasmania and southern New Zealand provide evidence for a rapid early or mid-20th century sea-level rise (Gehrels et al., 2008(Gehrels et al., , 2012. While potentially consistent with enhanced Northern Hemisphere land-based ice loss, questions remain over the existence and extent of this rise and the causal processes, in part due to perceived discrepancies with tide-gauge data from the region (Fadil et al., 2013). This paper therefore sets out to address the scarcity of Southern Hemisphere proxy records and establish whether the rapid rates of 20th century sea-level rise can be replicated and accounted for using the latest microfossil-based RSL reconstruction approaches. To this end, we present a new late Holocene RSL record from Mokomoko Inlet, near Bluff, southern New Zealand (Fig. 1), and enhance the temporal and vertical resolution of a previously published record from the nearby site of Pounawea (Gehrels et al., 2008).
As well as being few in number, existing proxy-based late Holocene RSL records from the south-west Pacific rely on limited quantification of the contemporary relationship between microfossil distributions and elevation above sea level (Gehrels et al., 2008(Gehrels et al., , 2012Figueira, 2012;Grenfell et al., 2012;Hayward et al., 2012). A comprehensive understanding of the elevation-dependent zonation of microfossils, e.g. foraminifera in intertidal environments, underpins quantitative proxy RSL reconstructions (Scott and Medioli, 1978;Birks, 1995;Hayward et al., 1999;Gehrels, 2000;Kemp and Telford, 2015;Williams et al., 2021). In recognition of the need for the development of more comprehensive training sets of modern samples, we enhance an existing foraminiferal dataset from Pounawea (Southall et al., 2006) with the addition of new samples from Mokomoko Inlet. In doing so, we seek to investigate the distribution of tidal-marsh foraminifera both locally and regionally, identify the controls on their distributions, and develop new transfer function models to underpin sea-level reconstructions from this important region. By applying knowledge of the contemporary distribution of foraminiferal assemblages and through Bayesian age-depth modelling and Gaussian process regression modelling, we aim to reconstruct RSL from fossil foraminifera preserved in tidal-marsh sediments at Mokomoko Inlet and Pounawea. By increasing the number and enhancing the resolution of RSL reconstructions in this understudied region, we endeavour to quantify rates of 19th and 20th century sealevel change and contribute to understanding the processes driving Southern Hemisphere sea-level rise.

Study sites
Mokomoko Inlet is situated on the southern side of the New River Estuary, approximately 14 km south-west of Invercargill, near the southern tip of New Zealand's South Island (Fig. 1). The inlet lies~110 km west of the previously studied site of Pounawea (Southall et al., 2006;Gehrels et al., 2008). The tidal range at the mouth of the New River Estuary ranges from 1.8 to 2.4 m during neap and spring tides, respectively, with highest astronomical tide (HAT) 1.5 m above mean sea level (MSL) (Todd, 2007;LINZ, 2019). Flows into and out of the estuary are dominated by tides, with the discharge rates through the 1-km-wide entrance channel of~1700-2300 m 3 s −1 being 30-40 times greater than fluvial inputs into the estuary (Todd, 2007). The embayment of Mokomoko Inlet, to the south of the estuary entrance channel, contains salt marshes and tidal flats. The inlet features a narrow bedrockframed mouth, <200 m wide, and a broader basin fringed by tidal marsh. Sedimentation in the centre of the inlet predominantly consists of sand transported from outside the estuary by tidal currents, with freshwater inputs limited to a few small streams and drainage channels (Thoms, 1981). Anthropogenic reclamation of the New River Estuary commenced in the mid-to late 19th century; however, these modifications have focused on the northern reaches of the estuary, with limited impacts on Mokomoko Inlet (NRETAC, 1973). The small tidal range, thin sedimentary sequence overlying incompressible basement rocks, limited freshwater input, stable geomorphic setting and the low likelihood of past changes in tidal range make the site suitable for sea-level reconstructions. This area of southern New Zealand is considered tectonically stable over multimillennial timescales (Beavan and Litchfield, 2012); however, there is uncertainty over the amount and spatial variability of tectonically driven vertical land motion over shorter timescales .  (Figueira, 2012;Figueira and Hayward, 2014), and Pounawea (Southall et al., 2006;Gehrels et al., 2008). Aerial photo of Mokomoko Inlet in panel c is from Bing Maps (https://www.bing.com/maps, image copyright DigitalGlobe, 2020). Panels d and e show the topography and vegetation zonation across the sampled area of marsh. See text for description of vegetation zones. Tape measure and flags in panel d mark the modern sample transect shown in e. Highest astronomical tide (HAT) and mean higher high water (MHHW) are indicated in panel e. [Color figure can be viewed at wileyonlinelibrary.com] The tidal marshes fringing Mokomoko Inlet display a clear vertical zonation in vegetation (Fig. 1d). The upper marsh (vegetation zone 3 in Fig. 1d) is characterized by Phormium tenax (New Zealand flax), Apodasmia similis (oioi or jointed wire rush) and Plagianthus divaricatus (saltmarsh ribbonwood). The mid-marsh (vegetation zone 2 in Fig. 1d) is composed of A. similis and Pl. divaricatus, while the lower marsh (vegetation zone 1 in Fig. 1d) contains only A. similis. Terrestrial vegetation above the limit of tides includes Leptospermum scoparium (manuka). The site lacks morphological evidence for erosion and reworking of sediment, with a gradual transition from low marsh to mudflat rather than an eroding cliff (cf. Figueira and Hayward, 2014).
The tidal marsh at Pounawea and its distribution of foraminifera in contemporary surface samples are described in detail by Southall et al. (2006). Further details concerning the stratigraphic record at Pounawea are provided by Gehrels et al. (2008).

Modern sample collection and preparation
To assess the distribution of tidal-marsh foraminifera with respect to elevation at Mokomoko Inlet, we collected 31 surface samples along a 45-m-long transect from the upper tidal flat to the zone of terrestrial vegetation ( Fig. 1; −46.537°, 168.282°). In keeping with other similar studies (e.g. Gehrels, 2000;Edwards et al., 2004;Williams et al., 2021), we focus on the upper half of the intertidal zone as we seek to develop a training set suitable for high-resolution sea-level reconstructions from organic salt-marsh sediments. Elevation is a surrogate variable, exerting no direct influence on species distributions. Rather, elevation is typically correlated with a number of other variables, including the frequency and duration of tidal flooding, vegetation type, salinity, grain size and organic carbon content (Wright et al., 2011). The transect covered a total elevation gradient of~1.1 m, with an average vertical spacing between samples of 0.035 m (range 0-0.08 m). The lack of geodetic benchmarks or tidal observations from the site resulted in differential GPS-derived sample elevations being expressed relative to a temporary benchmark. We convert these elevations to heights with respect to a tidal datum by identifying the highest occurrence of foraminifera, which approximates highest astronomical tide (HAT) (Gehrels et al., 2001;Charman et al., 2010;Wright et al., 2011). Further details on the use of tidal datums are provided in the section below.
Each 10-cm 3 sample consisted of the uppermost 1 cm of sediment, following a widely used sampling strategy designed to average out seasonal fluctuations in assemblages (Scott and Medioli, 1978;Horton and Edwards, 2003). This sampling strategy may miss the contribution of infaunal species that live at greater depths (Goldstein and Harben, 1993;Hayward et al., 2014;Chen et al., 2020), but ensures that significant changes in sea level have not occurred during the accumulation of the sediment. We added a Rose Bengal-ethanol solution to each sample within 24 h of collection to differentiate between living and dead foraminifera (Walton, 1952) and classed those with staining around the aperture and in the final chamber as living (Figueira, 2012). All samples were sieved, with the fraction between 63 and 500 μm retained and refrigerated in 5% ethanol before analysis. We used a wet splitter (Scott and Hermelin, 1993) to subdivide each sample into eight equal aliquots. Where possible, we counted a minimum of 100 dead specimens per sample under water. Where this number was not present in a 1/8 aliquot, we examined further splits until the desired minimum was reached or the entire sample was counted. We calculated densities based on count sizes and the number of aliquots examined. Some samples contained low densities and we only include modern samples with total dead counts >30 in our statistical analyses. While samples with low total counts may miss rare taxa, it is the abundant taxa that drive RSL reconstructions derived using transfer functions (see 'Statistical analysis of assemblages'). Consequently, reconstruction performance stabilizes at lower count sizes than are typically used in such studies (Payne and Mitchell, 2009;Kemp et al., 2020). Furthermore, our target count sizes are substantially larger than our accepted minimum (in contrast to the uniformly low count sizes explored by Kemp et al., 2020), providing a pragmatic approach that maximizes both accuracy and the quantity of samples included.
To compare the distribution of foraminifera at Mokomoko Inlet and other sites, and to develop a regional dataset, we incorporate previously published assemblage data from Pounawea (Southall et al., 2006) (Fig. 1). While Figueira and Hayward (2014) report the distribution of foraminifera in three transects from Waikawa Harbour, located between Mokomoko Inlet and Pounawea ( Fig. 1), we do not incorporate these samples due to substantial erosion and redistribution of sediments and foraminifera at this site, as noted by the original authors. At Pounawea, Southall et al. (2006) report one transect of 31 samples from 0 to 1 cm depth. Hayward et al. (2007) report nine further samples from lower elevations; however, these are from tidal flats rather than marsh environments. The calcareous species in these samples are not found in core samples and very low total counts of dead agglutinated specimens preclude their use here.
Species identifications follow Hayward and Hollis (1994) with updates to nomenclature following Hayward et al. (2020). While we have differentiated between Trochamminita salsa and Trochamminita irregularis in the Mokomoko Inlet surface samples following Callard et al. (2011), we group these species as Trochamminita spp. to enable comparison with the Pounawea modern and fossil datasets that do not make this subdivision. Our T. irregularis and Trochamminita spp. may also contain the newly recognized Pseudotrochamminita malcolmi (King, 2021).

Tidal data
We use the elevation of the highest occurrence of foraminifera as a common datum between sites (Wright et al., 2011). The highest occurrence generally lies close to HAT (Gehrels et al., 2001;Charman et al., 2010). As the tidal range differs between sites, we standardize elevations using HAT predictions from the TPXO8-Atlas global tidal model (Egbert and Erofeeva, 2010). To evaluate the performance of the model in this region, we first compare model predictions with data from the two closest tide gauges that have high-resolution data available from the Intergovernmental Oceanographic Commission Sea Level Station Monitoring Facility (www.iocsealevelmonitoring.org/). Puysegur Welcome Bay and Dunedin lie~140 and~190 km from Mokomoko Inlet, respectively (Fig. 1a). The tidal model successfully replicates data from these tide gauges, with decimetre-scale disparities and no need for changes in model phase or amplitude scaling (Supplementary Information S1). The predicted MSL to HAT range for Mokomoko Inlet, 1.49 m, also closely matches the 1.5-m range stated by Todd (2007) for the New River Estuary mouth. Consequently, we conclude that the tidal model can be applied to the field locations. The model predicts an MSL to HAT range of 1.19 m at Pounawea. To allow comparison of species' distributions between sites, we use these values to convert sample elevations to a Standardized Water Level Index (SWLI), where 100 SWLI is the highest occurrence of foraminifera and 0 SWLI is MSL.

Tidal-marsh stratigraphy and fossil foraminifera
Through the layers of sediment that accumulate over time and the microfossil assemblages contained within them, tidal marshes preserve a record of past changes in sea level (Gehrels et al., 2005;Barlow et al., 2013). We examined the stratigraphy of the tidal marsh at Mokomoko Inlet by cleaning the wall of a pre-existing~25-m-long trench, recently dug by a farmer to drain an adjacent field, in December 2006. We recorded the stratigraphy at nine equally spaced intervals using the Troels-Smith sediment classification scheme (Troels-Smith, 1955). Representative wide-bore gouge cores (10 cm in diameter) taken immediately adjacent to the logged sections provided material for subsequent analyses.
We subsampled the principal core, MKT-18.5, for foraminiferal analyses at 1-cm resolution, following the same preparation method as detailed above for the modern samples. We also counted foraminifera in 20 samples from six additional cores. We base our reconstructions on fossil samples with total dead counts exceeding 75 to ensure that the reconstruction precision is in line with the anticipated decimetre-scale RSL variability (Kemp et al., 2020).
The stratigraphy and fossil foraminiferal assemblages from the tidal marsh at Pounawea are detailed by Gehrels et al. (2008). With the exception of splitting two species of the genus Trochamminita (see above), the present study is consistent in approach with this earlier work.

Statistical analysis of assemblages
To investigate clustering in the modern foraminiferal assemblage data, we apply the Partitioning Around Medoids (PAM) algorithm (Rousseeuw, 1987;Kaufman and Rousseeuw, 1990;Kemp et al., 2012) in MATLAB v.2019b. Silhouette widths quantify how well each sample is classified and range from −1 (incorrect classification) to 1 (correct classification). We find the average silhouette width for 2-20 clusters, selecting the optimum number of clusters that provides the highest average silhouette width. We also visualize the relationships between samples using non-metric multidimensional scaling (NMDS) in the R package vegan (Oksanen et al., 2020). We choose Bray-Curtis dissimilarity and passively project modern sample elevation. The link between foraminiferal assemblages and an environmental gradient can be further described through constrained ordination. We use Detrended Canonical Correspondence Analysis (DCCA;ter Braak, 1986) in the CANOCO software package v.4.54 (ter Braak and Šmilauer, 2002) to quantify the variance in the dead assemblage dataset explained by the elevation gradient.
To enable the foraminiferal datasets to be used to reconstruct RSL change from fossil assemblages, we develop transfer functions to model the relationship between modern foraminifera and elevation (Imbrie and Kipp, 1971;Birks, 1995;Kemp and Telford, 2015). We explore the performance of models based both on a local training set containing data from a single site and a regional training set combining sites. While local models may provide lower sample-specific errors, regional models may provide closer analogues for fossil assemblages because they cover a wider range of environmental conditions (Gehrels et al., 2001;Wilson and Lamb, 2012;Watcham et al., 2013;Hocking et al., 2017;Rush et al., 2021). We first use DCCA to determine whether datasets show linear or unimodal relationships between assemblages and elevation. Birks (1995) suggests a threshold of 2 standard deviation units, with higher values indicating unimodal species responses. Based on this test, we employ linear Partial Least Squares (PLS; Wold et al., 1984) or unimodal Weighted Averaging (WA;ter Braak, 1987) or Weighted Averaging Partial Least Squares (WAPLS; ter Braak and Juggins, 1993) transfer function models in the C2 software package v1.7.3 (Juggins, 2007). Bayesian transfer functions offer an alternative approach and alleviate the need to assume linear or unimodal species responses (Holden et al., 2008;Cahill et al., 2016). While Bayesian approaches are not currently widely used for sea-level reconstructions, we apply Cahill et al.'s (2016) BTF package.
We assess transfer function performance using the crossvalidated (bootstrapped, 1000 cycles) correlation between observed and predicted elevations (r 2 boot ) and the root mean squared error of prediction (RMSEP). Following Juggins and Birks (2012), we review all samples with bootstrapped residuals greater than 2.5 standard deviations from the mean and remove them if there is additional evidence to suspect they are 'outliers' that do not reflect the long-term taxa-environment relationship.
We use the final accepted transfer functions to reconstruct the elevation at which each fossil sample from Mokomoko Inlet or Pounawea was deposited. The process of calibration provides palaeo-marsh surface elevation (PMSE) estimates with sample-specific errors. We convert these to estimates of RSL relative to present mean sea level by taking the PMSE away from the field elevation of the sample.
To assess whether modern training sets provide appropriate ranges of analogues for fossil samples, we calculate minimum dissimilarity coefficients (MinDC) using the squared chord distance metric (Birks, 1995). Following established convention we categorize all samples with MinDC values less than the fifth percentile of dissimilarity values in the modern training set as having a good modern analogue, those less than the 20th percentile as having a close modern analogue and all others as having a poor modern analogue (Garrett et al., 2013;Watcham et al., 2013). We supplement dissimilarity coefficients by visualizing taxonomic agreement between the modern and fossil samples using NMDS.
We test the significance of transfer function predictions by comparing the variance in the fossil assemblages explained by a transfer function reconstruction with the variance explained by 999 transfer functions trained on normally distributed random environmental variables (Telford and Birks, 2011). We perform this test using randomTF in the PalaeoSig R package (Telford and Trachsel, 2019).

Chronology
To determine the age of each layer of sediment in core MKT-18.5 and the age of the basal samples in other cores from Mokomoko Inlet, we use radiocarbon, radiocaesium, stable lead isotopes and pollen chronohorizons, and Bayesian age-depth modelling. Where possible, we radiocarbon dated fragile horizontally bedded terrestrial plant macrofossils, removing any visible roots or rootlets. Where such macrofossils were not available, we dated beetle carapaces. In an attempt to reduce temporal uncertainties, we divided some samples into two, analysed them separately and averaged the resulting ages. We report dates as conventional radiocarbon ages ( 14 C years before present) or the fractionation-corrected fraction modern carbon value (F 14 C) and calibrate to years CE using the ShCal20 calibration curve (Hogg et al., 2020).
We use the short-lived radionuclide 137 Cs, stable lead isotope ratios and a pollen chronohorizon to further constrain the age of the Mokomoko Inlet sedimentary sequence. The activity of 137 Cs was determined in contiguous 1-cm samples to identify the depth of the peak fallout, with analysis being undertaken using an Ortec HPGe well-type coaxial lowbackground intrinsic germanium detector at the CORiF laboratory, Plymouth, UK. We use stable lead isotopes to identify the onset of atmospheric pollution in the late 19th century and other ages related to both anthropogenic pollution and the 1815 CE eruption of Mount Tambora (Gehrels et al., 2008(Gehrels et al., , 2012. Lead isotope data were acquired using a Thermo Fisher Scientific Neptune multi-collector ICP-MS at the National Oceanographic Centre, Southampton, UK. We analysed pollen in contiguous 1-cm samples to identify the depth of the first occurrence of pine pollen and the presence of ruderal species, changes that can be linked to European settlement in the 19th century (Gehrels et al., 2008). Sample preparation for pollen analysis followed van der Kaars (1991) and identifications follow the Australasian Pollen and Spore Atlas (APSA Members, 2007).
We develop an age-depth model for core MKT-18.5 using a Bayesian framework in the R package rbacon (Blaauw and Christen, 2011;Blaauw et al., 2021).

RSL reconstructions
To reconstruct changes in RSL over time, we calibrate fossil foraminiferal assemblages from Mokomoko Inlet core MKT-18.5 and the previously published core, PW1, from Pounawea. Gehrels et al. (2008) report foraminifera from ã 40-cm-long core from the high marsh at Pounawea. In the original publication, the authors derived sea-level index points from nine precisely dated levels within the core. Here, we additionally calibrate foraminiferal assemblages from intervening levels, interpolating the corresponding ages from a new age-depth model based on the chronological data in the original publication.
We estimate rates of RSL change over time at Mokomoko Inlet and Pounawea using a Gaussian process regression model that considers both vertical and chronological errors assuming that they are both normally distributed, with temporal correlations in vertical uncertainties being approximated using an autoregressive process of the order 1 (Cahill et al., 2016;Gehrels et al., 2020). The temporal covariance is represented by a Matern function and all Gaussian process priors are estimated using a hyperparameter optimization scheme that is built in the MATLAB function fitgpr. A Monte-Carlo framework with 5000 iterations is used to estimate sea level change with its respective uncertainties at each site. The resulting 5000-member ensemble is used to assess statistics of rate changes in RSL.
To enable comparison between the two sites and with tide-gauge data, we correct sea-level records for changes in RSL resulting from glacial isostatic adjustment (GIA). Rather than using predictions from a single GIA model, these corrections use an ensemble of 5000 estimates, allowing for quantification of the GIA uncertainty . To understand the causes of 20th century sealevel changes, we also compare the proxy reconstructions and tide-gauge observations with sea-level budgets consisting of barystatic sea-level estimates from Frederikse et al. (2020), a sterodynamic estimate from the SODAsi.3 ocean reanalysis model (Giese et al., 2016)  Modern foraminifera: distribution and transfer function development

Distribution of foraminifera at Mokomoko inlet
The 31 surface samples from Mokomoko Inlet contain six species of foraminifera ( Fig. 2; Supplementary Information S2). Test concentrations are predominantly between 150 and 1000 per 10 cm 3 , with seven samples containing lower concentrations and insufficient dead counts (<30) to be included in our statistical analyses. Excluding these samples, the average count size is 123 and the maximum 181. The uppermost sample contains no foraminifera and we did not encounter any calcareous specimens in any sample. The six species display clear zonation across the marsh, with the PAM algorithm indicating the presence of three largely non-overlapping clusters in the dead assemblage data (Fig. 2). Cluster 1, containing six samples from the high marsh, is dominated by Trochamminita spp., which together make up 65-100% of the dead assemblage. Entzia macrescens and Trochammina inflata also reach their maximum abundances in this zone (11% and 4%, respectively). Haplophragmoides wilberti is the dominant mid marsh species, contributing between 50 and 80% of the dead assemblage in the 13 samples of cluster 2. In the low marsh, cluster 3 contains Miliammina fusca at abundances of 50-100%. Samples from the unvegetated mudflat contain a monospecific M. fusca assemblage, but with low concentrations and low total counts.
Stained foraminifera (assumed living at the time of collection) provide one-third of the total count on average (range: 0-67%). The living assemblage contains the same six species, with their distributions closely corresponding to the distribution of dead specimens (Fig. 2). Both E. macrescens and T. inflata are encountered at higher percentages in the living assemblage (up to 19 and 17%, respectively) than in the dead assemblage.
Elevation-dependent zonation in the regional dataset Foraminifera show vertical zonation in the Mokomoko Inlet surface samples, as they do in the dataset from Pounawea (Southall et al., 2006). Figure 3 compiles the modern assemblage data from these sites to summarize the regional distribution of the most abundant species with respect to elevation (SWLI). Of the 31 samples from Pounawea (Southall et al., 2006), 29 provide dead counts exceeding 30 specimens (average 147, maximum 253). The dead assemblage contains nine species (six agglutinated), with six (five agglutinated) exceeding 5% of the total dead count in at least one sample ( fig. 3 in Southall et al., 2006). The highest elevations are dominated by Trochamminita salsa, which probably also includes Trochamminita irregularis and is here referred to as Trochamminita spp. Samples from the mid-marsh at Pounawea feature abundant T. inflata, which reaches a maximum of 87%. By contrast, the species is rare at Mokomoko Inlet, reaching just 4%. Entzia macrescens is most abundant (40-60%) at the seaward end of the Pounawea transect, with M. fusca present at low percentages (average~5%) throughout the surface samples. There are no samples from lower elevations equivalent to those with high proportions of M. fusca that we report from Mokomoko Inlet.
The PAM algorithm separates the 53 samples of the regional training set into four clusters (Fig. 3). The samples contained within each cluster largely reflect the distribution of the dominant species, with cluster 1 containing the samples with abundant Trochamminita spp., cluster 3 containing the samples with dominant H. wilberti and so on. While there is some overlap between clusters, the influence of elevation is apparent. Samples in clusters 1 and 2 are restricted to SWLI values in excess of 90 and 73, respectively, while those in cluster 4 are strictly from elevations of <65 SWLI. NMDS results further underline the importance of elevation through the close alignment of the passively projected elevation gradient with the first axis (Fig. 3c). Constrained ordination quantifies the relationship with elevation; DCCA indicates that the elevation gradient explains 48.5% of the variance in the assemblage data from Mokomoko Inlet and 22.9% of the variance in the combined Mokomoko-Pounawea dataset.
While PAM clusters 1, 3 and 4 include samples from more than one site, cluster 2 contains samples from Pounawea only. These samples are dominated by T. inflata, a species infrequently encountered at Mokomoko Inlet. The separation of these samples is also apparent in the NMDS results, with samples from Pounawea clustering around T. inflata (Fig. 3c). The unusually low abundance of this species at Mokomoko Inlet may reflect the hypohaline (<30 psu) rather than euhaline (30-35 psu) nature of this site .

Transfer function development
We develop transfer functions based on both the local Mokomoko Inlet training set and the combined Mokomoko Inlet -Pounawea regional training set. DCCA confirms a unimodal relationship between assemblages and elevation for both training sets, supporting the use of WA and WAPLS models (Table 1; Birks, 1995). While Southall et al. (2006) reported a local training set from Pounawea, the removal of samples with low total counts leaves a total elevation range of just 0.27 m. The resulting local Pounawea transfer function (PLS, component 1) has little predictive power (bootstrapped r 2 = 0.23), emphasizing the need for a larger training set covering a greater elevation range.
The local Mokomoko Inlet WA and WAPLS transfer function display strong relationships between observed and predicted elevations, with no outliers with residuals exceeding 2.5σ ( Fig. 4a; Table 1). The total number of samples (24) is below most guidelines for the minimum required (~30-40) (e.g. Reavie and Juggins, 2011;Juggins and Birks, 2012;Watcham et al., 2013); however, these suggestions are predominantly based on diatom transfer functions with far greater species diversities. The low species diversity, combined with the acceptable performance and the close correspondence between species coefficients for dominant taxa in the local and regional models ( Supplementary Information S4), suggests that the local training set models the species-environment relationship effectively. We select the two-component WAPLS model as this provides a decrease in RMSEP exceeding 5% and an increase in r 2 boot . Due to the close similarity between the modern distributions of T. salsa and T. irregularis (Fig. 2), grouping these two species as Trochamminita spp. does not negatively impact on the performance of the local Mokomoko Inlet transfer function and indeed results in a slightly lower RMSEP ( Supplementary Information S3).
The combined Mokomoko Inlet and Pounawea regional model also provides a strong relationship between observed and predicted elevations, but with larger maximum residuals, including one sample in the chosen two-component WAPLS model exceeding 2.5σ from the mean (Fig. 4b). This Pounawea sample lies on a local topographic high (see sample 15 in fig. 3 of Southall et al., 2006) and we choose to remove it due to the anomalously high abundance of M. fusca, the highest in any sample at the site. We note that the authors of the original study also considered this and three other samples as outliers. We do not remove the other three samples as they do not exceed our chosen thresholds for consideration. At Mokomoko Inlet, such high M. fusca abundances are only found at elevations below 73 SWLI, 20 SWLI units below the Pounawea sample. With this sample removed, we select the twocomponent model based on the >5% decrease in RMSEP over the one-component model and enhanced performance over the WA models (Table 1). Unusually, the regional RMSEP is lower than the local, suggesting no loss of precision with the larger training set.
We note that the E. macrescens species coefficient in the regional training set is low compared with the observed distribution of the species (Supplementary Information S4). This reflects rising percentages at the lowest end of the Pounawea transect and the lack of lower samples from this site in the training set (Fig. 3). The discrepancy between the species coefficients in the local and regional datasets is not a cause for concern as E. macrescens is rare in fossil sequences (Gehrels et al., 2008) and uncertainty over the species' optimum elevation is therefore unlikely to influence RSL reconstructions.
A Bayesian transfer function trained on the 53-sample regional training set did not perform as well as the WAPLS model, with a larger RMSEP and residual structure indicating overprediction at lower elevations and underprediction at higher elevations ( Supplementary Information S5). Consequently, we choose to use the more widely applied WAPLS approach for fossil reconstructions.
The accepted transfer function models, the local Mokomoko Inlet and regional (Mokomoko plus Pounawea) WAPLS models, display slightly larger RMSEP values than the local Pounawea model of Southall et al. (2006). This apparent decrease in performance reflects the larger sampling range covered by our training sets, which is likely to more faithfully represent species' distributions. When expressed as a percentage of the sampled elevation gradient, both of our transfer functions outperform the earlier model (~8% for the models presented here, 20% for the Pounawea model).

Stratigraphy
The cleaned wall of the trench, ranging from 0.2 m in height at the landward end to 0.6 m at the seaward end,  Figure 4. Performance of a. the local Mokomoko Inlet and b. the regional transfer function models. The left-hand panels show observed elevation against transfer functionpredicted elevation; the right-hand panels show residuals (predicted minus observed). WAPLS: weighted averaging partial least squares; r 2 boot : bootstrapped r 2 ; RMSEP: root mean square error of prediction; SWLI: standardized water level index. [Color figure can be viewed at wileyonlinelibrary.com] provided a cross-section of the marsh stratigraphy ( Fig. 5a; Supplementary Information S6). The base of the trench is marked by a layer of subangular greywacke boulders, cobbles and very coarse gravel over lithified greywacke. At the seaward end of the trench, the basal unit is overlain by silty sand and organic sandy silt, while at the landward end, the lowest fine-grained sediments are brown muddy peat or dark brown peaty mud. The latter is found along the length of the trench, overlying the organic sandy silt at the seaward end. The surficial sediments consist of greybrown peaty mud. Together, the two peaty mud units arẽ 0.3 m thick, with both containing abundant roots including the living roots of the current salt-marsh vegetation. None of the contacts between successive fine-grained sedimentary layers are abrupt or show evidence for erosion, and we did not encounter any coarser-grained layers that might indicate deposition by storms or other extreme waves. Core MKT-18.5, recovered from immediately adjacent to the trench, is 0.51 m long and consists of organic sandy silt overlain by 0.34 m of peaty mud (Fig. 6b).

Geochronology
Nineteen radiocarbon age determinations, including duplicate ages on two samples, provide constraints on the ages of the sediments at Mokomoko Inlet (Table 2). Radiocarbon ages on basal samples range from 859 ± 22 14 C years to modern, with no clear relationship between age and elevation (Fig. 5a). The nine samples from core MKT-18.5 also show age reversals, with four modern ages likely to result from contamination from younger root material which could not always be conclusively removed.
To address the issues encountered with the radiocarbon results, we applied other complementary dating approaches. Low unsupported 210 Pb concentrations preclude the use of this chronometer, as previously identified at Pounawea (Gehrels et al., 2008). 137 Cs concentrations, however, reach a clear peak at 0.065 m depth (Fig. 5b), which we attribute to 1965 CE, allowing 2 years for interhemispheric mixing and full atmospheric homogenization in the Southern Hemisphere (Hancock et al., 2011;Gehrels et al., 2012;Foucher et al., 2021). Changes in the 206 Pb/ 207 Pb curve are less pronounced than at Pounawea (Gehrels et al., 2008) or in Tasmania (Gehrels  Table 2. b. Chronohorizons in the master core, MKT-18.5. c. Lead isotope data from core MKT-18.5 in the context of source material including andesites and basalts from the Egmont and Taupo volcanic zones on the North Island (Graham et al., 1992;Price et al., 1992;Huang et al., 2000), Tambora lava (Turner and Foden, 2001), pre-contamination peat from Lynch's crater (Kylander et al., 2007), ores from Broken Hill and Mount Isa (Chiaradia et al., 1997;Vallelonga et al., 2002), and 1990s atmospheric values from Hamilton, Auckland and Christchurch (Bollhöfer and Rosman, 2000). d. Bayesian age-depth model for core    (Fig. 5b). A three-isotope plot reveals an excursion at 0.415 m towards values known from the 1815 CE Tambora eruption (Fig. 5c), as also recorded at Pounawea (Gehrels et al., 2008).
Indicators of disturbance (both bracken spores and charcoal particles) indicate the base of core MKT-18.5 postdates Polynesian settlement; however, the lack of diagnostic European pollen indicators implies a pre-European settlement age (Supplementary Information S7). While pollen is sparse between 49 and 29 cm, assemblages from 29 to 27 cm are probably pre-European. Introduced pine pollen is first encountered at 26-27 cm, alongside ruderal disturbance indicators. We assign a calendar age of 1880 ± 10 CE at this depth, following historical records of initial European settlement of the Southland Plains (Supplementary Information S7).

Biostratigraphy
Samples from core MKT-18.5 contain six species of foraminifera (Fig. 6), all of which are also found in the surface samples from the site. Total dead counts exceed 75 in every sample and average 124. The basal grey-brown organic sandy silt contains a mixed assemblage of M. fusca, T. salsa and H. wilberti, which each exceed 50% of the assemblage in different samples from this unit. Above 0.38 m, T. salsa dominates, exceeding 70% of the dead assemblage in all samples between 0.38 and 0.09 cm. The uppermost 0.1 m sees an increase in H. wilberti, which exceeds 60% of the assemblage at the top of the core, and M. fusca, which reaches~10%. Test densities average around 90 per cubic centimetre (ml −1 ), with peaks exceeding 300 ml −1 and lower densities in the uppermost and lowermost 0.07 m.
Calibration of the dead assemblages in core MKT-18.5 using the local and regional transfer functions provides palaeomarsh surface elevation predictions that are around mean higher high water (MHHW) at the base of the core, rising to near HAT for the T. salsa-dominated interval above 0.38 m (Fig. 6). The increase in H. wilberti and M. fusca at the top of the core drives a decline in palaeomarsh surface elevation of~0.1 m. The predictions from the local and regional models are in strong agreement, with an average difference of 0.01 m. Uncertainties (2σ) associated with the local model average 0.16 m, while those of the regional model average 0.14 m.
The local Mokomoko Inlet training set fails to provide good modern analogues for 32 of the 50 fossil samples in core MKT-18.5, a figure that is reduced to 14 samples with the use of the regional training set (Fig. 6). These 14 samples are all located below 0.34 m depth, with five below 0.43 m lacking even close analogues in the regional training set. Joint NMDS ordination of the modern and fossil samples supports the similarity in assemblages for much of the core, with only the base of the core plotting away from the scatter of modern samples (Supplementary Information S8). The samples without good modern analogues have mixed assemblages of species that are found at different elevations in the modern environment, raising the possibility that they have been influenced by post-depositional processes. In particular, the abundance of the high marsh T. salsa despite the inorganic sedimentation and the co-dominance of the low marsh or mudflat-inhabiting M. fusca could suggest alteration of the assemblage through infaunality or differential preservation of tests.
Our palaeomarsh surface elevation reconstruction explains 60% of the variance in the fossil assemblages in core MKT-18.5, but narrowly fails to outperform the suite of 999 random transfer functions (95th percentile = 62%). Such significance test failures have previously been reported in sea-level reconstructions from high-marsh deposits because of the  limited variability in reconstructed elevations as the marsh surface keeps pace with sea-level rise (Kemp et al., 2013). Consequently, we do not view the significance test result as any cause for concern. Basal samples from   Fig. 5) yielded low total counts (two, three and nine specimens, respectively), preventing their use for reconstructing sea level. Cores MKT-5 and MKT-12 also had low counts in their basal samples; however, in both cases, overlying samples provide some corroboration of the environment in which the basal sedimentary units were deposited (Fig. 7). In both instances, the lowermost samples contain T. salsa, which is also dominant in the overlying samples (>70%), indicating deposition in a high marsh environment. The low densities in the basal samples may indicate proximity to the highest occurrence of foraminifera; therefore, rather than using the transfer function predictions, we assign an indicative meaning of HAT with an uncertainty of 0.15 m. Miliammina fusca dominates the basal sandy silt in core MKT-15.5; however, the overlying sample contains a mixed assemblage with a higher proportion of T. salsa. The disparity in assemblages, the presence of T. salsa in organic-poor sediments, the substantial difference in predicted marsh surface elevation, the low density for both samples and the poor analogue for the upper sample suggest possible postdepositional modification through infaunality or differential test preservation. Consequently, we do not use this core for reconstructing sea level.

Sea-level reconstructions Mokomoko Inlet
We reconstruct RSL at Mokomoko Inlet by combining the regional PMSE predictions (subtracted from the field elevation of each sample) with the age-depth model or the individual basal radiocarbon ages. We favour the regional over the local model due to the larger proportion of fossil samples that have good modern analogues in this training set and the greater precision. We limit the reconstruction to the organic section of core MKT-18.5 and the organic-rich basal samples that show the least evidence for post-depositional alteration (see 'Biostratigraphy' above). We correct each sea-level point for GIA at a rate of −0.04 mm a −1 (95% interval −0.21 to 0.26 mm a −1 )  and use the Gaussian process regression model to quantify rates of change over time for the period after 1850 CE.
Our new reconstruction reveals that GIA-corrected relative sea level (GCSL) was a few tens of centimetres below present in the 13th century CE and fell to~0.3-0.6 m below present by the late 15th to early 17th century CE (Fig. 7a). GCSL rose bỹ 0.4 m between 1850 CE and the end of the record in~2002 CE (Fig. 7b). Sea levels rose more gradually during the mid-to late 19th century, with rates averaging 1.6 mm a −1 (95% interval −1.1 to 3.9 mm a −1 ) between 1850 and 1900 CE. The rate of rise increased during the early 20th century, peaking in 1940 CE at a rate of 5.0 (3.4-9.5) mm a −1 , before declining during the latter half of the century. From 1980 CE until the end of the record, rates averaged 1.7 (−0.9 to 4.9) mm a −1 . Supplementary Information S9 provides the vertical position of sea level and the corresponding age range for each analysed sample.

Pounawea
To improve the RSL reconstruction from Pounawea, we recalibrate the fossil foraminiferal assemblages reported by Gehrels et al. (2008) using the regional transfer function and develop a new Bayesian age-depth model using chronological information from the original paper (Supplementary Information S10). Our reconstruction is consistent with the previously published sea-level curve, but calibrating assemblages at additional depths provides a more finely resolved record (Figs. 6 and 7; Supplementary Information S11). The reconstruction again fails the significance test due to the prevalence of high marsh samples (SWLI explains 53% of the variance vs. 62% explained by 999 random transfer functions); as discussed in 'Biostratigraphy' above, this is not a cause for concern.
Corrected for GIA at a rate of −0.04 (−0.21 to 0.26) mm a −1 , our improved reconstruction suggests RSL was stable a few tens of centimetres below present between the start of the record in the late 15th to mid-17th century CE (Fig. 7b). In close correspondence with the Mokomoko Inlet reconstruction, RSL rose gradually in the 19th century, before rapidly accelerating in the early 20th century and slowing during the late 20th century. Reconstructed 19th century rates are slower at Pounawea than at Mokomoko Inlet, averaging 0.03 (−1.8 to 1.8) mm a −1 , while the early 20th century acceleration was more rapid and the mid-century peak slightly higher at 6.0 (4.0-10.7) mm a −1 . The fastest reconstructed rates of sea-level rise occurred in 1946 CE. From 1980 until the end of the record, rates averaged 0.7 (−2.1 to 2.9) mm a −1 . The Pounawea and Mokomoko Inlet reconstructions provide a consistent picture of a rapid mid-20th century sea-level rise and a late 20th century decline in rates but differ slightly in the timing of the initiation of the 19th or early 20th century acceleration and the total late Holocene rise in sea level.

Reconciling proxy and instrumental records
Rapid sedimentation rates, mediated by rising RSL and increased accommodation space, facilitate a high temporal resolution for the late 19th and 20th centuries at both Mokomoko Inlet and Pounawea, overlapping with the period covered by the longest tide-gauge records. In Fig. 7(c) we compare the proxy reconstructions with the instrumental record from Dunedin (location in Fig. 1), which provides an annual record going back to 1899 CE (Hogarth, 2014;Frederikse et al., 2020). The tide gauge lies~190 km from Mokomoko Inlet and 90 km from Pounawea. While the Bluff tide gauge is situated within 10 km of Mokomoko Inlet, it only includes six annual mean values from before 1986 and we therefore focus comparisons on the more distant tide gauge. While Gehrels et al. (2008) excluded the Dunedin record from their comparison due to suspected wharf subsidence (Hannah, 2004), Denys et al. (2020) subsequently confirmed its stability. The proxy and instrumental records concur that RSL rose over the 20th century, but the Dunedin tide gauge records a smaller total rise and lacks clear evidence for more rapid rates in the early to mid-20th century. This discrepancy could indicate that the tidal-marsh records overestimate RSL rise or that the locations do not share a common RSL history because of varying driving processes. Fadil et al. (2013) suggest the Pounawea sequence could be influenced by sediment compaction, resulting in overestimation of RSL. However, compaction modelling demonstrates that thin tidal-marsh sequences overlying lithified rock or well-consolidated sediments, such as those at Mokomoko Inlet and Pounawea, are not significantly affected, with compaction likely to induce millimetre-scale postdepositional changes at most (Brain et al., 2012). King et al. (2021) raise the possibility of the small modern training set used by Gehrels et al. (2008) influencing the accuracy of the reconstruction. With the addition of new modern samples from Mokomoko Inlet and the development of a robust regional transfer function, we can also discount this possibility. The close correspondence in reconstructions between Mokomoko Inlet and Pounawea and our focus on organic sequences dominated by mid-to high-marsh sediments that closely track sea-level rise further argue against the proxy records overestimating RSL rise.
Changes in RSL not only reflect changes in sea-surface height, but also changes in land level. Differences in RSL histories between the tidal marsh sites and Dunedin could, therefore, reflect differential vertical land motion at the tide gauge and proxy locations. While uplift at Dunedin could reconcile the records, continuous global positioning system (cGPS) time series indicate subsidence at a rate averaging 1.12 ± 0.22 mm a −1 over the last 20 years (Denys et al., 2020). If this subsidence rate was constant it would fully account for the locally observed rates of 20th century RSL rise and it therefore appears likely that the short cGPS record overestimates the long-term trend (Denys et al., 2020;King et al., 2021). InSAR data suggest a high degree of spatial variability in subsidence around Dunedin with lower rates, −0.75 ± 0.14 mm a −1 , averaged across a 5-km area around the tide gauge and higher rates, −2.40 ± 1.40 mm a −1 , along the city's waterfront (Levy et al., 2020). A 5.5-year cGPS time series from Bluff, close to our Mokomoko Inlet site, suggests uplift of 0.4 mm a −1 , albeit with a large uncertainty of 0.7 mm a −1 (Beavan and Litchfield, 2012). Extrapolation of this shortterm pattern of vertical land motion again does not serve to resolve the difference with the instrumental record, but rather increases the disparity. Furthermore, New Zealand's southern coastline is generally considered tectonically stable over multimillennial timescales (Beavan and Litchfield, 2012;Clement et al., 2016;Ryan et al., 2021), suggesting that recent uplift may be a transient feature.
While stable over long timescales, Mokomoko Inlet, Pounawea and Dunedin are potentially affected by coseismic and postseismic vertical land motion associated with major earthquakes (magnitude 7 + ) generated by subduction along the Puysegur Trench to the south-west of the South Island (Denys et al., 2020). With substantial distances separating the locations, differential vertical motions (both co-and postseismic) associated with major subduction earthquakes are feasible. The 2009 Dusky Sound earthquake (M w 7.8), with an epicentre 300 km west of Dunedin, 250 km from Pounawea and 160 km from Mokomoko Inlet, uplifted both Dunedin and Bluff by a few millimetres (Beavan et al., 2010). Land-level changes associated with earlier major earthquakes in 1938, 1939, 1945and 1979are unknown (Denys et al., 2020. Upper plate structures associated with the Otago fold-thrust belt provide a further potential source for vertical land motion at Dunedin and Pounawea (Langridge et al., 2016); however, the New Zealand GeoNet Quake Search (quakesearch.geonet. org.nz) indicates no earthquakes exceeding magnitude 5 have occurred in this region over the 20th century. Indeed, the Settlement Fault, which intersects the Pounawea estuary~2 km from the studied core, is not known to have ruptured in the last millennium (Hayward et al., 2007). Coseismic, postseismic and interseismic deformation remain largely unknown contributors to RSL changes in this region over centennial to multicentennial timescales.
While we have been unable to fully resolve the discrepancies in rates of sea-level rise between the proxy records and the tide-gauge observations, replication of the rapid 20th century rates at a second tidal marsh provides additional weight to support the findings of Gehrels et al. (2008) and counters suggestions that the RSL rise has been overestimated. Longer cGPS time series and ongoing tide-gauge measurements may in time help to resolve the differences and to identify the importance of transient vertical land motion, including the influence of regional earthquake deformation.

Mechanisms driving 20th century RSL change in southern New Zealand
We attempt to account for the pattern and rates of GIA-corrected RSL rise in the Mokomoko Inlet and Pounawea records by deriving a sea-level budget for each location. These budgets combine barystatic sea-level estimates (Frederikse et al., 2020), sterodynamic estimates (Giese et al., 2016) and the inverse barometer effect (Slivinski et al., 2019). The budgets for Mokomoko Inlet and Pounawea show a high degree of similarity to our GIA-corrected RSL reconstructions, including the higher rates centred on the 1940s and lower rates in the early and late 20th century (Fig. 8c). The largest differences between the proxy records and the budgets are seen between 1900 and 1920, with the budgets lying above the upper margin of the Gaussian process 95% intervals, but still well within the vertical uncertainty of the individual sealevel index points. Individual budget components also carry greater uncertainties towards the start of the 20th century (Giese et al., 2016;Frederikse et al., 2020). The agreement between the proxy reconstructions and the budgets suggests that barystatic and sterodynamic processes are sufficient to explain the rapid rates of 20th century sea-level change reconstructed for southern New Zealand.
Globally, enhanced rates of 20th century sea-level rise are well documented, with particularly rapid rates (~2.4 mm a −1 ) centred on the 1930s and 1940s reflecting an above-average barystatic contribution, specifically enhanced melt from glaciers and the Greenland Ice Sheet (Marzeion et  Marzeion, 2021). Barystatic estimates are, however, marginally lower in southern New Zealand than the global average during this period. Despite Northern Hemisphere ice losses generally resulting in faster rates of rise in far-field locations due to the spatially varying fingerprint of ice melt, this was counteracted in New Zealand by the gravitational fingerprints of local ice loss (Frederikse et al., 2020). Nevertheless, barystatic rates at our sites were still above the long-term average during the early to mid-20th century (Fig. 8a).
The contribution of Antarctic melt to barystatic sea-level change before the satellite era remains almost entirely unknown (Galassi and Spada, 2017;Rignot et al., 2019). Observations suggest some mass loss in the mid-to late 20th century from West Antarctica (Smith et al., 2017), with more rapid rates during 21st century (Galassi and Spada, 2017;Shen et al., 2018;Shepherd et al., 2018). Our budget follows Frederikse et al. (2020) in assuming a small and unchanging contribution from Antarctica of 0.05 ± 0.04 mm a −1 before 1993, based on the compilation of Adhikari et al. (2018). The rapid rates of sea-level rise that we have reconstructed from New Zealand are unlikely to reflect enhanced Antarctic ice melt. If Antarctic melt had substantially contributed to these rates, far-field Northern Hemisphere sites would be expected to show more rapid rises and such rises have not been observed (Gehrels and Woodworth, 2013;Gehrels et al., 2020). Nevertheless, quantifying Antarctic contributions to sea-level rise over the pre-satellite era remains of great importance for contextualizing current enhanced ice melt and should be the focus of future research.
In addition to enhanced Northern Hemisphere ice loss, the early to mid-20th century also coincides with a period of rapid rates of sterodynamic sea-level rise around New Zealand, reflecting both increases in upper ocean heat content and atmospheric and oceanic circulation changes (Giese et al., 2016;Zanna et al., 2019). Indeed, it is the sterodynamic component that chiefly accounts for the particularly rapid rates of change in our budgets for Mokomoko Inlet and Pounawea (Fig. 8a). Mapping this component reveals a region with extremely rapid rates of rise centred on southern New Zealand between 1940 and 1950 (Fig. 8b). A range of interconnected atmospheric and oceanic phenomena may account for this pattern, including the Interdecadal Pacific Oscillation (IPO) and El Niño Southern Oscillation (ENSO), alongside changes in the position and strength of Southern Hemispheric westerlies (Swart and Fyfe, 2012;Dangendorf et al., 2019) and possible changes in the strength of the subtropical gyre (Goring and Bell, 1999;de Lange, 2001;Salinger et al., 2001;Albrecht et al., 2019). The most rapid rates of sea-level rise in our proxy reconstructions coincide with the transition from a positive to a negative phase of the IPO, which was characterized by an increase in sea-surface temperatures and consequently sea levels around New Zealand (Goring and Bell, 1999;de Lange, 2001;Hannah and Bell, 2012). The IPO also modulates the strength and frequency of ENSO events (termed El Niño and La Niña), with resulting changes in the frequency of storm surges and changes in winds, currents and wave climate providing further mechanisms for sea-level fluctuations. While La Niña events did occur in the 1950s, potentially enhancing sea-level rise, the 1930s and 1940s were marked by weaker ENSO activity (Goring and Bell, 1999) and the contribution of this climatic phenomenon to the mid-20th century rise in southern New Zealand remains somewhat uncertain. Nevertheless, persistent El Niño conditions, coinciding with lower sea-surface temperatures associated with a positive IPO phase from 1976 to the late 1990s, are likely to have contributed to a late 20th century deceleration in sea-level rise (Bell et al., 2000), as seen in our proxy reconstructions.
Our proxy records from Mokomoko Inlet and Pounawea both end in the early 2000s; consequently, neither record currently documents 21st century sea-level trends. More rapid sea-level rise might be anticipated during the first few decades of this century due to the switch to a negative IPO phase and enhanced barystatic contributions related to anthropogenically forced ice melt (Frederikse et al., 2020).

Conclusions
In this paper, we explore the elevation-dependent distribution of foraminifera in modern tidal marsh sediments from southernmost New Zealand as a basis for reconstructing late Holocene RSL change. We develop a regional transfer function that provides crucial underpinning for enhanced sea-level reconstructions from the Southern Hemisphere, both in this paper and potentially through future studies. Such reconstructions are critical for quantifying the contributions of the different processes driving sea-level change in the Common Era.
Mokomoko Inlet, a tidal marsh close to Bluff, provides an ideal location for late Holocene and historical RSL reconstructions. We achieve these through calibrating fossil foraminiferal assemblages using the regional transfer function, comprehensive dating and Bayesian age-depth modelling. GIA-corrected sea level was~0.3 m below present in the 13th century CE and fell to~0.5 m below present by the late 15th to early 17th century CE. Sea levels rose gradually during the mid-to late 19th century, before rapidly accelerating during the 20th century, reaching a peak rate of 5.0 (95% interval 3.4-9.8) mm a −1 in 1940 CE. Slower rates characterize the late 20th century, averaging 1.7 (−0.9-5.1) mm a −1 from 1980 to the end of the record in the early 2000s.
We also revisit an existing record from Pounawea (Gehrels et al., 2008), applying the new regional transfer function, enhancing the temporal resolution. and applying newer age, GIA and Gaussian process regression modelling techniques. Together, the Pounawea and Mokomoko Inlet reconstructions provide independent yet consistent histories of late Holocene RSL change in southern New Zealand. Although these proxy reconstructions cannot yet be fully reconciled with observations from the closest long tide gauge at Dunedin, spatial differences in co-and postseismic vertical land motions could provide feasible explanations. Furthermore, the sea-level reconstructions from Mokomoko can be comprehensively explained through comparison with current best estimates of contributors to sea-level change over the 20th century. Rapid rates of RSL rise in southern New Zealand in the early to mid-20th century are consistent with the coincidence of enhanced melt from glaciers and the Greenland Ice Sheet and particularly fast sterodynamic sea-level rise associated with atmospheric and oceanic circulation changes. The late 20th century deceleration is similarly consistent with persistent El Niño conditions and lower sea-surface temperatures associated with a positive IPO phase.

Supporting information
Additional supporting information can be found in the online version of this article.This article includes online-only Supplemental Data. Figure S1: Comparison of tide-gauge data and tidal model predictions for a. Puysegur Welcome Bay and b. Dunedin. Table S2: Percentage abundance of foraminifera (unstained and stained) in the Mokomoko Inlet surface samples. We provide separate counts for Trochamminita salsa and T.
irregularis; these species are lumped as Trochamminita spp. in the local and regional transfer functions. Supplementary Information S3 compares the performance of local transfer function models with Trochamminita lumped or split. Figure S3: Comparison of the performance of two transfer functions based on the local Mokomoko Inlet training set with Trochamminita spp. either a. lumped or b. split. The left-hand panels show observed elevation against transfer function predicted elevation; the right-hand panels show residuals (predicted minus observed). Figure S4: Predictions of species' coefficients (elevation optima and tolerances) in the local Mokomoko Inlet and regional transfer function models. We group Trochamminita salsa and T. irregularis as T. spp. SWLI: Standardised Water Level Index; WAPLS C2: Weighted Averaging Partial Least Squares component 2. Figure S5: Performance of a regional Bayesian transfer function developed using the BTF package (Cahill et al., 2016). a. Observed elevation against Bayesian transfer function predicted elevation. We reject this model due to the larger RMSEP and residual structure indicating overprediction at lower elevations and underprediction at higher elevations. b. Species response curves. Figure S6: Location and photographs of the trench originally dug by a farmer through the marsh at Mokomoko Inlet to drain an adjacent field. The stratigraphy in Figure 5 is based on the cleaned wall of this trench. In c, Roland Gehrels recovers a core from adjacent to the trench. Circles in d note the locations of particular cores, with along-trench distances indicated. Figure S7: Summary of pollen analyses from core MKT-18.5 from Mokomoko Inlet. Sample preparation for pollen analysis followed van der Kaars (1991) and identifications follow the Australasian Pollen and Spore Atlas (APSA members, 2007). Calculation of pollen percentages followed standard New Zealand palynology practice (e.g. Newnham et al., 1995) Figure S8: Simultaneous NMDS ordination of the regional training set and the fossil assemblages from Mokomoko Inlet, with elevation (SWLI) passively projected. Modern samples indicated by dots; fossil samples are joined by a blue line showing the stratigraphic order, with the core top at the thin end of the line and the bottom at the thick end. We perform NMDS using the R package vegan (Oksanen et al., 2019), with Bray-Curtis dissimilarity. Table S9: Sea-level index points for Mokomoko Inlet, without corrections for glacial isostatic adjustment. MinDC is the Minimum Dissimilarity Coefficient calculated using the squared-chord distance metric. We provide MinDCs for both the local Mokomoko Inlet training set and the regional training set. We use the 5th percentile of the modern training set dissimilarities as the boundary between good and close analogues (2.84 for the local training set, 7.52 for the regional training set) and the 20th percentile as the boundary between close and poor modern analogues (13.99 for the local training set, 30.46 for the regional training set). Close modern analogues are italicised, poor modern analogues are in bold and italicised. We provide the RSL data in the HOLSEA workbook format at DOI: 10.6084/m9.figshare.19134335. Figure S10: Age-depth model for Pounawea core PW1 developed using rBacon (Blaauw and Christen, 2011;Blaauw et al., 2021) and the SHCal20 radiocarbon calibration curve (Hogg et al., 2020). Table S11: Sea-level index points for the Pounawea core of Gehrels et al. (2008) without corrections for glacial isostatic adjustment. MinDC is the Minimum Dissimilarity Coefficient calculated using the squared-chord distance metric. We provide MinDCs for both the regional training set presented in this paper and, for comparison, the PW06 training set of Southall et al. (2006). Two close (rather than good) modern analogues are italicised.