Characterisation and quantification of organic carbon burial using a multiproxy approach in saltmarshes from Aotearoa New Zealand

Blue carbon ecosystems, such as saltmarshes, play a vital role in mitigating climate change by sequestering atmospheric carbon dioxide and storing it as buried organic carbon for centuries to millennia. While there are international methodologies for generating blue carbon credits through coastal wetland restoration, their application in Aotearoa New Zealand is limited due to insufficient data on saltmarsh carbon stocks, accumulation rates and the processes governing long-term carbon preservation. To quantify these metrics, we examined 45 sediment cores collected from five saltmarsh sites in Aotearoa New Zealand. The cores were analysed for elemental composition, stable isotopes, and lipid biomarkers. These data were collected using a range of techniques, including X-ray fluorescence (XRF), Ramped-Pyrolysis Oxidation-Accelerator Mass Spectrometry (RPO-AMS), and Pyrolysis-Gas Chromatography-Mass Spectrometry (Py-GC-MS). Results show high variability in soil organic matter properties, carbon stocks (41.3±9.4 to 92.3±66.2 Mg C ha⁻¹; mean ± SE), and accumulation rates (0.46±0.02 to 1.53±0.09 Mg C ha⁻¹ yr⁻¹; mean ± SE). Stable isotope and lipid biomarker results indicate substantial contributions from saltmarsh vegetation to the organic carbon pool. Results suggest that plant-derived organic carbon is preserved in the oldest basal sediments. Our findings highlight that spatial variability must be considered when conducting carbon assessments in saltmarsh ecosystems. Further research is required to determine the environmental drivers that influence long-term carbon storage and to improve the accuracy of blue carbon assessments in Aotearoa New Zealand.

1 Introduction

Coastal wetlands, such as saltmarshes, mangroves, and seagrass meadows, sequester atmospheric carbon dioxide (CO₂) and store it as buried organic carbon (OC; “blue carbon”) over centuries to millennia (Chmura et al., 2003; Mcleod et al., 2011). These ecosystems can accumulate up to 53.65 Tg OC yr⁻¹ and account for 30 % of the total carbon burial in ocean sediments (Wang et al., 2021). This storage potential has led to high public and private interest in protecting and restoring these ecosystems for their climate mitigation potential (e.g., Ministry for the Environment, 2022, 2024; Ross et al., 2024), and blue carbon credit methodologies have been developed for the voluntary carbon market (Friess et al., 2022; Lovelock et al., 2023a; Needelman et al., 2018). As a result, research that aims to quantify the carbon sequestration potential of blue carbon ecosystems (BCEs) has increased over the past decade (Howard et al., 2023; Macreadie et al., 2019). Despite these efforts, major gaps in the data required to fully characterise these systems remain, especially in temperate regions of the Southern Hemisphere (Bertram et al., 2021; Howard et al., 2023; Macreadie et al., 2019, 2021).

Saltmarshes occur at the interface between terrestrial, marine, and estuarine settings and accumulate OC that is: (a) produced in-situ by saltmarsh plants (autochthonous sources) and (b) derived from terrestrial and marine organisms that live outside the saltmarsh and are transported and deposited at the marsh surface by riverine runoff and tidal inundation (allochthonous sources; Howard et al., 2014; Middelburg et al., 1997). These OC sources often mix with siliciclastic minerals that are also transported via fluvial and coastal currents and are deposited to form minerogenic soils (Howard et al., 2014; Middelburg et al., 1997; Saintilan et al., 2013). Given the complex combination of autochthonous and allochthonous sediment contributions to the below-ground carbon pool, improving our understanding of these sources and their preservation mechanisms is critical for accurately estimating sequestration potential and informing blue carbon credit methodologies.

The average soil carbon accumulation rate (CAR) for saltmarshes in Aotearoa New Zealand (NZ) is 0.89 Mg C ha⁻¹ yr⁻¹ (Bulmer et al., 2024), which is significantly smaller than the global mean estimates that range between 1.67 and 2.45 Mg C ha yr⁻¹ (Chmura et al., 2003; Ouyang and Lee, 2014; Wang et al., 2021). Attempts to characterise the source of OC and its preservation characteristics in BCEs in NZ are also limited (e.g., Bulmer et al., 2020; Pérez et al., 2017; Sikes et al., 2009; Thomson et al., 2025). Here we examine five coastal saltmarsh systems at three locations that span 6.56° of latitude from Rangaunu Harbour in Northland (34°58^′ S, 173°13^′ E) to Pāuatahanui Wildlife Reserve in Wellington (41°06^′ S, 174°54^′ E). Our aim is to better characterise these saltmarsh systems via three primary research objectives: (i) quantify the carbon stocks and accumulation rates across a range of saltmarsh habitats, (ii) assess the source and preservation characteristics of buried organic material, and (iii) synthesise results to guide future research.

2 Study sites

Five saltmarsh sites in three geographic locations in the North Island of NZ were selected to capture a range of geomorphic and environmental settings.

2.1 Okatakata Islands, Omaia Island and Awanui, Rangaunu Harbour, Northland

Omaia and Okatakata Islands are situated within the 9700 ha Rangaunu Harbour in Northland, which features tidal flats colonised by mangroves and saltmarshes, and three main channels: Kaimaumau, Awanui, and Pukewhau (Fig. 1; Heath et al., 1983). Stopbanks around the southern shores and tributary streams were first built in 1916 and prevent spring tides and tidal surges from extending inland (Cathcart, 2005). Okatakata consists of two islands, forming a 38 ha saltmarsh. Omaia is a 50 ha drained saltmarsh used as pastureland since 1937, with drains and stopbanks installed to prevent tidal flooding (Land Information New Zealand, 2018). Isolated patches of saltmarsh (Awanui) are present immediately south of Omaia. Historical imagery shows that Awanui and Okatakata have remained undisturbed by human activities since at least the early 1940s.

2.2 Robert Findlay Wildlife Reserve, Pūkorokoro-Miranda, Waikato

Robert Findlay Wildlife Reserve (Robert Findlay) spans approximately 27 ha on the Firth of Thames' western coastline, within the 8500 ha Pūkorokoro-Miranda coastal wetland (Fig. 1), which is designated a wetland of “international significance” (Gerbeaux, 2003). The wetland, located on a chenier plain, was used for farming since 1865, and a floodgate and drains were installed in the early 1900s when stilt ponds were created for shell extraction (Queen Elizabeth II National Trust, 1992; Woodley, 2016). By 1980, the drainage canals and floodgates were no longer maintained and gradually filled with silt, enabling the saltmarsh to regenerate. A conservation covenant was established in 1988 (Queen Elizabeth II National Trust, 1992).

2.3 Pāuatahanui Wildlife Reserve, Wellington

Pāuatahanui Wildlife Reserve (Pāuatahanui) is situated within the Pāuatahanui inlet, which forms the northern branch of the Te Awarua-o-Porirua Harbour, located 30 km north of central Wellington, and has extensive areas of saltmarsh approximately 50 ha in size (Fig. 1). The saltmarsh formed after a M_w 8.2 earthquake on the Wairarapa fault in 1855, which caused widespread uplift around the Wellington region (Grapes and Downes, 1997, 2010). Historical records indicate that prior to this earthquake, the site of the present-day saltmarsh comprised subtidal or tidal flat environments (McManaway and Gaz, 1852; Park, 1841). The marsh likely began to form immediately following the uplift event and was well established by ca. 1865 (Stephenson, 1986). Drainage canals were built across sections of the saltmarsh and the area was used as pastureland for cattle and sheep until the 1980s (Sheehan, 1988). In 1980, NZ Forest and Bird replanted four hectares of saltmarsh, and in 1984, the NZ Government established the Wildlife Management Reserve and began restoration on the remaining 46 ha (Conwell, 2010; Guardians of Pāuatahanui Inlet, 2021).

Figure 1 Study sites and sample locations within (a) Rangaunu Harbour, where Oma, Oka and Awa represent study sites Omaia, Okatakata and Awanui, respectively; (b) Robert Findlay, where core locations are denoted as Puk; and (c) Pāuatahanui, where core locations are represented as Pau. Imagery is sourced from the LINZ Data Service and licensed for reuse under the CC BY 4.0 licence.

3 Methods

3.1 Field sampling

Prior to sample collection, ground-based vegetation surveys were used to classify low, mid and high marsh zones for all sites based on saltmarsh plant community composition (after King, 2022). Historical aerial imagery was also examined to delineate former marsh zones across Omaia. Regions with higher former vegetation density were characterised as high and mid marsh, and lower density zones as low marsh. Detailed documentation of vegetation was done in 50×50 cm quadrats prior to sediment coring. Dominant vegetation types are summarised in Table S1 (Supplement).

Twelve sediment cores were collected from Pāuatahanui saltmarsh in November and December 2021. During fieldwork in January 2022, 33 sediment cores were collected from Omaia (9), Okatakata (9), Awanui (6), and Robert Findlay (9) (Fig. 1 and Table S1). All cores were collected with a gouge auger (6 cm diameter; 50 cm length), which recovers a cylindrical sediment core with minimal compaction (Smeaton et al., 2020). Cores were placed in PVC half-pipes with ice packs and transported to Earth Sciences New Zealand (ESNZ), Lower Hutt, NZ, and stored in a refrigeration facility at 4 °C. Cores were described following the Troels-Smith (1955) sediment classification system. One core from each marsh zone from Omaia, Okatakata and Robert Findlay and two cores from each marsh zone from Pāuatahanui were sampled in 2 cm depth increments from the top of the core down to 50 cm, and after that in 5 cm intervals down to the base of the core (after Howard et al., 2014). The remaining cores were sub-sampled at 5 cm intervals between 0 and 20 cm and in 10 cm intervals from 20 cm to the base. A total of 97, 25, 56, 75, and 61 samples were collected from Pāuatahanui, Robert Findlay, Awanui, Okatakata and Omaia cores, respectively. In cases where the base of the saltmarsh deposit was difficult to identify, a microscope was used to determine whether saltmarsh foraminifera were present or not (after King et al., 2024).

3.2 Elemental and stable isotope analysis of Total Organic Carbon (TOC) and Total Nitrogen (TN)

First, large roots and aboveground biomass were manually picked from surface samples to avoid biasing average soil organic matter (OM) properties. We note that samples were not size fractioned as this study focuses on bulk soil OC, which includes belowground living plant biomass (e.g., small rootlets and rhizomes; Macreadie et al., 2017). All samples were then weighed, freeze-dried and weighed again, and homogenised using a ball mill.

Three hundred and fourteen samples were analysed for total organic carbon (TOC; wt %) and total nitrogen (TN; wt %) concentrations and their stable isotope composition (δ¹³C_org and δ¹⁵N) at the Stable Isotope Laboratory at ESNZ. TOC and δ¹³C_org were determined on acidified samples (treated with 10 % HCl for 12 h) by elemental analysis isotope ratio mass spectrometry (EA-IRMS) using an EA Eurovector 3000 and Elementar Isoprime model. TN and δ¹⁵N values were analysed on unacidified samples (Smeaton et al., 2024; Sollins et al., 1999). Internal reference standards for δ¹³C_org (Cane Sugar −10.3 ‰, beet sugar −24.6 ‰ and EDTA −31.1 ‰) and δ¹⁵N (Leucine 2.0 ‰, EDTA 0.58 ‰ and Caffeine −7.8 ‰) were run every 10 samples. δ¹³C_org and δ¹⁵N values are reported in permil (‰) relative to the Vienna Pee Dee Belemnite (VPDB) standard and AIR, respectively. The C:N ratio is reported as the molar ratio of TOC to TN. The analytical precision of the measurements is ±0.2 wt % for TOC, ±0.1 wt % for TN, ±0.2 ‰ for δ¹³C and ±0.3 ‰ for δ¹⁵N.

Dry bulk density (DBD; g cm⁻³) was calculated by dividing the mass of the dry sample by the sample volume following standard methodologies (Howard et al., 2014). Organic carbon density (CD; g C cm⁻³) was calculated by multiplying bulk density by the OC content for each sample depth increment (wt %; Howard et al., 2014). TOC stocks (Mg C ha⁻¹) for each core were calculated by integrating the depth intervals (2, 5 or 10 cm) over the depth range of the core.

3.3 Chronology

We use a published age model from core PauM1, which was previously collected from Pāuatahanui (King et al., 2024), to interpolate the age of sediments in core Pau-HM3 and estimate CARs. The age model for PauM1 was developed using a Bayesian framework in the R package rplum (Blaauw et al., 2024) and assumes the base of the Pāuatahanui saltmarsh formed in 1855 following uplift during a M_w 8.2 earthquake on the Wairarapa fault (Grapes and Downes, 1997; Stephenson, 1986).

Lead isotope data were used to produce new age-depth models for three cores: Okatakata (Oka-MM1), Awanui (Awa-MM2) and Robert Findlay (Puk-MM1). Each core was sub-sampled in 1 cm increments, and the number of samples varied from 9 to 16 based on the core length and stratigraphy. Lead isotope data were generated using gamma and alpha spectrometry conducted at the radio-isotope facility at the New Zealand Institute of Public Health and Forensic Science, Christchurch, NZ. All samples were measured with gamma spectrometry, and every second sample in each core was measured with alpha spectrometry. For gamma spectrometry, sediment samples were packed into petri dishes and left to equilibrate for three weeks and then analysed to detect radionuclide activity to include ²¹⁰Pb, ¹³⁷Cs, ²²⁸Ra and ²²⁶Ra (Arias-Ortiz et al., 2018; Goldstein and Stirling, 2003). For alpha spectrometry, the samples were first processed to prepare the granddaughter ²¹⁰Po source, and the activity of ²¹⁰Po was then measured to calculate excess ²¹⁰Pb activities. Decay of excess ²¹⁰Pb activity (half-life = 22.5 years; Appleby, 1997; Arias-Ortiz et al., 2018) and ¹³⁷Cs discharge peak in ∼1965 (Goff and Chagué-Goff, 1999) were used to determine the rate of sediment accumulation. Age-depth models were generated using rplum (Blaauw et al., 2024).

CARs for each core were calculated by dividing the TOC stock for each depth interval by the corresponding age as per the rplum age-depth models. CARs are presented for each age-depth interval with the 95 % confidence interval (CI). Mean CARs for each marsh core were then calculated and are presented as mean ± SE.

¹⁴C dating was attempted on sieved sedge and rush fragments (>1 mm) from basal samples from cores Oka-MM1, Awa-MM2 and Oma-MM3 at the Rafter Radiocarbon Laboratory at ESNZ. However, the reported calibrated ages for all basal samples were modern (post-bomb) and unsuitable for inclusion in the age-depth models and for estimating a date of original saltmarsh establishment at Omaia (Albot et al., 2026).

3.4 X-ray fluorescence (XRF)

Elemental abundances in 314 samples were measured using an Olympus Vanta M-series XRF portable scanner at Victoria University of Wellington, Wellington, NZ. The scanner has a 9 mm-diameter primary beam connected to a workstation that allows remote operation. Approximately 2–3 g of dried and homogenised sediment was placed in plastic tubes (5 mL; 15 mm diameter), resulting in sediment thickness in the tube of >10 mm. Each sample was measured using the standard “Geochem 3-Beam” method (50, 40 and 10 kV beams set for 30, 30 and 40 s “live time” respectively) built into the scanner, which is designed to detect and quantify major and some trace elements, including Al, Ca, Fe, K, Si, Mn, S, Sr, Ti, Zn, Cu, Pb, Zr, Nb and Rb. Measurements were taken at the start and end of each analytical period for nine external standards provided by the United States Geological Survey (USGS; AGV-2 Andesite; BHVO-2 Hawaiian Basalt; COQ-1 Carbonatite; W-2 Diabse; SGR-1 Green River Shale; SCo-1 Cody Shale) and the Geological Survey of Japan (GSJ; JR-2 Igneous; JG-2 Igneous; JF-2 Igneous). Results are reported as absolute concentrations in parts per million (ppm). The measurements were used to calculate elemental ratios relevant to OM sources and preservation. Specifically, Ca, Sr and S as marine sources (e.g., carbonates and sulfates), with Ca and Sr also reflecting biogenic carbonate; S as indicative of reducing conditions/anoxic environments that result from intrusion of sulphate-rich marine waters into organic sediments where sulphate-reducing bacteria oxidise OM; Ti, Al, Fe, Si, and K as terrestrial/lithogenic indicators from weathering of continental silicate rocks; Fe and Mn as markers of redox processes, for example, sulphate-reduction intensifying biogeochemical cycling of metals; and Zr:Rb as a grain size proxy (coarse-clay ratio), as Zr is predominantly found in coarser sediments and Rb in clays (Croudace and Rothwell, 2015; Ewers Lewis et al., 2019; Kelleway et al., 2017; Naeher et al., 2013).

3.5 Organic carbon fingerprinting

3.5.1 Lipid extraction and Gas Chromatography-Mass Spectrometry (GC-MS) analysis

Selected samples from Pāuatahanui (Pau-HM1 n=5, Pau-MM4 n=6, Pau-LM4 n=6), Okatakata (Oka-HM2 n=5, Oka-MM1 n=5, Oka-LM1 n=6) and Awanui (Awa-MM2 n=9) were analysed to determine their lipid biomarker compositions. Samples were selected along core profiles at intervals where pronounced shifts or distinct changes were observed in TOC contents, C:N ratios, δ¹³C_org and δ¹⁵N trends. Analyses were carried out in the ESNZ/VUW Organic Geochemistry Laboratory at ESNZ following the methodology described in Naeher et al. (2012, 2014) with some modifications, as reported in Verret et al. (2025).

The GC-MS data, focused on n-alkanes and steroids, were interpreted with Agilent MassHunter/Chemstation software based on relative retention times and diagnostic mass spectra. Several indices based on n-alkane distributions, such as the carbon preference index (CPI), odd-over-even predominance (OEP), the average chain length ratio (ACL), and the relative contribution of aquatic plants relative to terrestrial biomass (P_aq index), as well as C₂₈ and C₂₉ stanol-sterol ratios, were calculated as described below.

CPI provides an estimate for the predominance of odd-numbered over even-numbered carbon chains and was calculated for C₂₃-C₃₃ carbon homologues using Eq. (1) as reported in Wang et al. (2003):

$$\begin{array}{}\text{(1)}& \text{CPI}={\displaystyle \frac{\sum \text{oddC}n}{\sum \text{evenC}n}}\end{array}$$

CPI has been commonly used as a proxy to estimate the degree of OM degradation or determine dominant OM sources. Odd-numbered n-alkane chains dominate in fresh biomass and recent sediments, and diagenetic alteration of OM results in preferential decay of odd-chain-length n-alkanes (Bray and Evans, 1961; Meyers and Ishiwatari, 1993). High CPI values indicate fresh OM, whereas values <1 indicate a high degree of microbial degradation or thermal maturation (Cranwell, 1981; Eglinton and Hamilton, 1967). Values close to 1–2 indicate partial OM degradation by microorganisms (Jaffé et al., 2001; Tanner et al., 2010; Zhao et al., 2024). Emergent and submerged/floating aquatic plants typically exhibit CPI values >3 (Bray and Evans, 1961; Eglinton and Hamilton, 1967; Jiménez-Morillo et al., 2024).

Similar to CPI, OEP also represents the predominance of odd-numbered over even-numbered carbon chains. This index is calculated using Eq. (2) from Zech et al. (2010).

$$\begin{array}{}\text{(2)}& \text{OEP}={\displaystyle \frac{{\mathrm{C}}_{\mathrm{27}}+{\mathrm{C}}_{\mathrm{29}}+{\mathrm{C}}_{\mathrm{31}}+{\mathrm{C}}_{\mathrm{33}}}{{\mathrm{C}}_{\mathrm{28}}+{\mathrm{C}}_{\mathrm{30}}+{\mathrm{C}}_{\mathrm{32}}}}\end{array}$$

High OEP values >1 have been interpreted to represent fresh, undegraded OM or terrestrial plant sources, while lower OEP values are indicative of OM degradation or less terrestrial inputs (Cranwell, 1981; Eglinton and Hamilton, 1967; Wang et al., 2003, 2015; Zech et al., 2010).

ACL is the weighted average of carbon chain lengths for the long-chain n-alkanes detected in the C₂₇-C₃₁ range and was calculated following Eq. (3) of Poynter and Eglinton (1990):

$$\begin{array}{}\text{(3)}& \text{ACL}={\displaystyle \frac{\sum [\left(\mathrm{27}\times {\mathrm{C}}_{\mathrm{27}}\right)+\left(\mathrm{29}\times {\mathrm{C}}_{\mathrm{29}}\right)+\left(\mathrm{31}\times {\mathrm{C}}_{\mathrm{31}}\right)]}{\sum ({\mathrm{C}}_{\mathrm{27}}+{\mathrm{C}}_{\mathrm{29}}+{\mathrm{C}}_{\mathrm{31}})}}\end{array}$$

Variations in ACL can represent changes in vegetation types and have also been partly attributed to changes in prevailing temperature and/or moisture of the surrounding environment (Derrien et al., 2017; Poynter et al., 1989; Zhou et al., 2010). For example, the predominance of C₂₇ and C₂₉ n-alkanes is characteristic of rush, sedge, shrub, and tree species, while C₃₁ and C₃₃ are more abundant in grasses and herbs (Eley et al., 2016; Zech et al., 2010).

P(aqueous), referred to as P_aq, is calculated using Eq. (4) from Ficken et al. (2000):

$$\begin{array}{}\text{(4)}& P_{\mathrm{aq}}={\displaystyle \frac{{\mathrm{C}}_{\mathrm{23}}+{\mathrm{C}}_{\mathrm{25}}}{{\mathrm{C}}_{\mathrm{23}}+{\mathrm{C}}_{\mathrm{25}}+{\mathrm{C}}_{\mathrm{29}}+{\mathrm{C}}_{\mathrm{31}}}}\end{array}$$

Higher P_aq values indicate a higher proportion of submerged and emergent aquatic vascular plants (macrophytes) and wetter conditions (Ficken et al., 2000). Low P_aq values of ≤0.25 indicate high contributions of higher terrestrial vascular plants, whereas high P_aq values >0.4 indicate dominant contributions of marine and freshwater macrophytes (Ficken et al., 2000; Zhao et al., 2024). Sikes et al. (2009) further differentiate between emergent macrophytes, typically ranging between 0.4 and 0.6, and submerged and floating macrophytes exhibiting values >0.6.

Finally, the stanol-sterol ratio was calculated following Eq. (5) from Naafs et al. (2019):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \text{stanol}-\text{sterol ratio}\\ & ={\displaystyle \frac{{\mathrm{C}}_{\mathrm{28}}\text{stanol}+{\mathrm{C}}_{\mathrm{29}}\text{stanol}}{{\mathrm{C}}_{\mathrm{28}}\text{stanol}+{\mathrm{C}}_{\mathrm{29}}\text{stanol}+{\mathrm{C}}_{\mathrm{28}}\text{sterol}+{\mathrm{C}}_{\mathrm{29}}\text{sterol}}}\end{array}\end{array}$$

C₂₈ and C₂₉ sterols are derived mainly from higher plants (Gaskell and Eglinton, 1976; Volkman, 1986). In anaerobic and reducing conditions, sterols are reduced to stanols due to microbial activity. This makes the stanol-sterol ratio a good indicator of diagenetic degradation of plant OM, where higher stanol-sterol ratios indicate more degraded material (Naafs et al., 2019; Wakeham, 1989).

3.5.2 Ramped-Pyrolysis Oxidation-Accelerator Mass Spectrometry (RPO-AMS)

RPO-AMS analysis was run at the Rafter Radiocarbon Laboratory at ESNZ, following the methodology described in Ginnane et al. (2024), with minor modifications (i.e., sieving to obtain a homogenous, non-biased sample). This analysis thermochemically separates OC in a particulate sample to create a radiocarbon profile of its compositional constituents (Ginnane et al., 2024; Rosenheim et al., 2008). Analyses were performed on one basal sample each from Okatakata (Oka-MM1, 28–30 cm depth) and Awanui (Awa-MM2, 90–95 cm depth). Δ¹⁴C values and conventional radiocarbon ages (CRAs) are reported as defined by Stuiver and Polach (1977) and fraction modern (F_m) values are reported as defined by Donahue et al. (1990).

The Δ¹⁴C data for each oxidised pyrolytic split is used to calculate the relative proportion of syndepositional (i.e., modern) OC versus recalcitrant carbon (i.e., older/reworked carbon) based on the modified isotopic mixing model by Broz et al. (2023) provided in Eq. (6), below.

$$\begin{array}{}\text{(6)}& C_{\mathrm{modern}}=\text{TOC}\left({\displaystyle \frac{\mathrm{\Delta }{}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{split}}-\mathrm{\Delta }{}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{last}\mathrm{\_}\mathrm{split}}}{\mathrm{\Delta }{}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{modern}}-\mathrm{\Delta }{}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{last}\mathrm{\_}\mathrm{split}}}}\right)\end{array}$$

where C_modern is the modelled fraction of syndepositional carbon, TOC is the total OC content (wt %) of the bulk sample, Δ¹⁴C_split is the measured Δ¹⁴C value of each pyrolytic split, Δ¹⁴C_modern is a typical value for a modern post-bomb OC endmember (where Δ¹⁴C is assumed to be 0 ‰ in materials deposited in the last 2 kyr), and Δ¹⁴C_{last_split} is the final Δ¹⁴C pyrolytic split value where all labile carbon components are considered to be degraded based on the Pyrolysis-Gas Chromatography-Mass Spectrometry (Py-GC-MS) composition results as described below.

3.5.3 Pyrolysis-Gas Chromatography-Mass Spectrometry (Py-GC-MS)

Py-GC-MS analyses imitate the ramped pyrolysis process without oxidation to determine the molecular distributions that correspond to the radiocarbon measurements, which provide information about the origin, relative quantities, and degradation states of different OM sources (Ginnane et al., 2024). Py-GC-MS analysis was undertaken in the ESNZ/VUW Organic Geochemistry Laboratory at ESNZ, following Maier et al. (2025). These analyses were performed on one basal sample each from Okatakata (Oka-MM1, 28–30 cm depth) and Awanui (Awa-MM2, 90–95 cm depth).

Pyrolytic compounds for each split were grouped into the following categories (after Carr et al., 2010; Kaal et al., 2020; Maier et al., 2025; Zhang et al., 2019; and references therein): n-alkanes, with < C₂₁ and ≥ C₂₁ representing marine and terrestrial vegetation sources, respectively; polysaccharide derivatives (e.g., furans, furaldehydes and related compounds) from plant pigments; thiophenes derived from sulphur compounds; phenols derived from terrestrial plant lignin; polycyclic aromatic hydrocarbons (PAHs) as indicators of terrestrial carbon sources; and nitrogen (N)-containing compounds (e.g., benzonitrile, indole and related compounds) as indicators of proteinaceous aquatic microorganisms or microbial biomass. The remaining compound classes, such as cyclic alkanes and alkylbenzenes, are considered undiagnostic because they are primarily derived from recalcitrant OM (e.g., Ginnane et al., 2024; Maier et al., 2025).

3.6 Statistical analysis

3.6.1 Soil organic matter properties

R Studio (version 4.2.1) was used to test whether variability in measured soil variables across five different study sites and between individual cores collected from three distinct vegetation zones within sites was statistically significant. First, Levene's test for equal variance and Shapiro-Wilk's test for normality were run on all datasets. Because the normality and equal variance assumptions were not met, non-parametric Kruskal-Wallis and post-hoc Dunn's tests were used for pairwise comparisons. To account for the increased likelihood of Type I errors (incorrectly finding a statistically significant difference when there isn't one) associated with multiple comparisons, the Benjamini-Hochberg correction was applied to control the false discovery rate.

3.6.2 Principal Component Analysis (PCA) and hierarchical clustering

Principal Components Analysis (PCA) and hierarchical clustering were conducted on two separate datasets that included elemental, stable isotope, and XRF data, with variations based on the presence or absence of lipid biomarker indices. PCA was performed in R Studio (version 4.2.1) using the FactoMineR package. Where results showed low inter-variable correlation, thereby limiting the effectiveness of PCA, variables with minimal association with PC1 and PC2 were excluded, leading to more interpretable PCs. Specifically, the quality of representation (Cos2), which identifies the contributions of the variables on each PC, and the loadings matrix, which indicates the strength and direction of the relationship (positive or negative), were used to interpret the strength and direction of relationships.

Hierarchical clustering utilising Ward's criterion was employed to complement PCA and provide an alternative perspective on data clustering. Analyses were performed in R Studio using the hclust function. Welch's t-tests were conducted to assess the statistical significance of differences between cluster means and the overall mean for each variable. Additionally, pairwise t tests with the Benjamini-Hochberg correction were performed to compare the means of variables between the identified clusters.

4 Results

4.1 Stratigraphy and sedimentology

Depth of refusal, typically indicating the base of saltmarsh sediments, ranged from a minimum thickness of 5 cm for low and high marsh cores at Robert Findlay to >95 cm for a mid marsh core at Awanui. Marsh sediments at Pāuatahanui comprised herbaceous peat or organic-rich silts in the top 5–15 cm, underlain by silts and clays (Fig. 2a). Shelly sands and clays, interpreted as pre-marsh sediments, were observed at refusal depths across all cores. Most marsh sediments at Robert Findlay generally comprised sandy/silty peat in the top 5 cm, underlain by shelly silts and sands in high marsh zones and shelly clay in low and mid marsh zones (Fig. 2b). At Omaia, cores consisted of a thin layer of topsoil underlain by layers of intermixed silty sands, silts and clays with occasional thin peaty layers (Fig. 2c). Sediments at Okatakata included a 10–15 cm thick upper unit of sandy/silty peat and organic-rich silts that sit on top of silty clays (Fig. 2d). Marsh sediments at Awanui typically comprised peaty/organic-rich silts and clays within the top 10 cm, underlain by organic-rich clays and silt (Fig. 2e, f). However, several cores collected from low marsh areas with juvenile vegetation (Awa-LM1 and Awa-LM2) were entirely composed of organic-rich sands. Most cores collected from low marsh areas at all sites were bioturbated. Examples of detailed stratigraphic logs following the Troels-Smith (1955) classification are provided for cores Pau-HM3, Puk-MM1, Oma-MM3, Oka-MM1, and Awa-MM2 in Figs. S1–S5 (Supplement).

Figure 2 Simplified stratigraphy for sediment cores from (a) Pāuatahanui, (b) Robert Findlay, (c) Omaia, (d) Okatakata, (e) Awanui, and (f) simplified stratigraphy and photograph of core Awa-MM2. The Troels-Smith (1955) classification is simplified to peaty sediment and minerogenic lithologies in this schematic diagram.

4.2 Bulk soil organic matter variables and carbon stocks

Statistically significant differences (p<0.05) in TOC, TN, CD and DBD were observed between sites and in at least one marsh zone at each site, although no consistent patterns emerged regarding which zones had higher values. Statistical results for Levene's, Shapiro-Wilk's and Kruskal-Wallis tests are presented in Table S2. Tables S3–S7 provide the results of the post-hoc Dunn tests (Supplement). Mean carbon stocks ranged from 41.3±9.4 Mg C ha⁻¹ at Robert Findlay to 92.3±66.2 Mg C ha⁻¹ at Awanui (Table 1; Fig. 3). In most marsh zones, mean carbon stocks exceed median values, demonstrating a positive skew arising from the sampling of locations with higher carbon stocks. Table S8 (Supplement) summarises the results by marsh zone for each site. Kruskal-Wallis results indicate no significant difference in carbon stocks across the five sites (N=45; p=0.48).

Table 1 Mean (± SE) values for Total Organic Carbon (TOC), Total Nitrogen (TN), Dry Bulk Density (DBD), Carbon Density (CD), and carbon stocks (mean ± SE) for each study site.

Figure 3 (a) Mean carbon stocks (Mg C ha⁻¹ ± SE) at each study site, and (b) mean carbon stocks per marsh zone for each study site. The black dotted line represents the median value.

4.3 Chronology

Total ²¹⁰Pb activities for all three age-depth models decline from the surface and reach supported levels (i.e., equilibrium with parent isotope; ≤10²¹⁰Pb Bq kg⁻¹) at 40 cm for Awanui and 20 cm at Okatakata, but do not reach supported levels at Robert Findlay. ¹³⁷Cs peaks at Robert Findlay, Okatakata and Awanui were evident between 9–10, 2–3, and 7–8 cm, respectively (Fig. 4b, d, f), and correlate to the ∼1965 peak fallout in NZ (Goff and Chagué-Goff, 1999). The mean 95 % CI is 26.2 years at Robert Findlay, 51.2 years at Okatakata and 24.9 years at Awanui. Robert Findlay shows age uncertainty <30 years for the top 10 cm of the core, which increases up to 47 years at the base of the marsh deposit (Fig. 4a). Age uncertainty at Okatakata is below 30 years for sediments deposited down to 5 cm and increases to 55 years at the marsh base (Fig. 4c). Age uncertainty in the Awanui core for sediments deposited between 0 and 20 cm is less than 20 years. This uncertainty increases to between 20–32 years below 20 cm (Fig. 4e).

Figure 4 Bayesian age-depth models for (a) Robert Findlay (core Puk-MM1), (c) Okatakata (core Oka-MM1), and (e) Awanui (core Awa-MM2). The y axis represents the depth of the cores (cm), and the x axis depicts the estimated age (years in Common Era; CE). The red dotted line represents the model mean, and the grey shaded area represents the 95 % confidence interval (CI). The dotted boxes and shading represent total ²¹⁰Pb (Bq kg⁻¹) activity as measured by alpha and gamma spectrometry. ¹³⁷Cs peak is plotted as a calendar date (1965±2). (b), (d) and (f) show ¹³⁷Cs levels for each core.

4.4 Carbon Accumulation Rates (CARs)

Soil CARs ranged from a minimum of 0.37 Mg C ha⁻¹ yr⁻¹, CI [0.35, 0.43], in basal sediments at Okatakata, to a maximum of 2.16 Mg C ha⁻¹ yr⁻¹, CI [1.96, 2.44], in surface sediments at Awanui (Fig. 5). Robert Findlay exhibited the highest mean CAR (1.53±0.09 Mg C ha⁻¹ yr⁻¹; mean ± SE), followed by Awanui (1.32±0.01 Mg C ha⁻¹ yr⁻¹), Pāuatahanui (0.64±0.02 Mg C ha⁻¹ yr⁻¹), and Okatakata (0.46±0.02 Mg C ha⁻¹ yr⁻¹).

Figure 5 Core soil CARs (Mg C ha⁻¹ yr⁻¹ ± mean 95 % CI) obtained from Bayesian age-depth models generated using rplum (Blaauw et al., 2024): Awanui core Awa-MM2; Okatakata core Oka-MM1; and Robert Findlay core Puk-MM1. Pāuatahanui core Pau-HM3 was interpolated based on the age-depth model in King et al. (2024).

4.5 Stable isotope analysis of carbon and nitrogen

δ¹³C_org, δ¹⁵N and C:N ranges for the study sites are provided in Table 2. δ¹³C_org versus C:N scatter plots (Fig. 6) show that the sample distributions across all sites fall predominantly within the C₃ plant range, freshwater dissolved OC (DOC) and particulate OC (POC) sources. Omaia and Okatakata show contributions from two additional sources of OC: freshwater algae and marine DOC. Omaia exhibits the widest range of carbon sources.

Figure 6 Scatterplots of δ¹³C_org (‰) and C:N ratios for (a) Pāuatahanui, (b) Robert Findlay, (c) Omaia, (d) Okatakata and (e) Awanui. Note that the x and y scales differ based on site-specific ranges. The coloured boxes represent the typical δ¹³C_org and C:N ranges of the different sources of organic inputs to the coastal environment – C₃ terrestrial plants, freshwater DOC, marine DOC, freshwater algae, marine algae, and freshwater POC. The ranges have been compiled from various studies and presented by Lamb et al. (2006).

Table 2 δ¹³C_org, δ¹⁵N and C:N ratio ranges at the five study sites.

4.6 X-ray fluorescence (XRF)

XRF analyses quantified the following major and trace elements in all cores: Mn, Al, Fe, Zr, Ca, Si, K, Sr, Nb, Rb, Ti, Zn, S, Pb, and Cu (Albot et al., 2026). PCA and hierarchical cluster analysis were used to examine the stratigraphy and relationships between XRF measurements and other datasets. These results are shown in Fig. 9 and discussed in Sect. 5.2 and 5.3.

4.7 Lipid biomarkers

Distribution of biomarkers, ratios and indices

The distribution of n-alkanes in the apolar biomarker fractions ranges from C₁₈ to C₃₃ and shows dominance of either mid-chain (mid-molecular weight; C₂₁-C₂₅) or long-chain (high molecular weight; C₂₆-C₃₃) n-alkanes across all sites. There is a higher relative abundance of odd-carbon n-alkanes than even-carbon n-alkanes (Fig. 7). Pāuatahanui has higher relative abundances of long-chain n-alkanes (65 %) than mid-chain n-alkanes (33 %). In comparison, at Okatakata and Awanui, the mid-chain n-alkanes (50 % and 75 %, respectively) exhibit a higher relative abundance than the long-chain n-alkanes (49.5 % and 25 %, respectively). Short-chain n-alkanes (< C₂₁) contributions were 2 % at Pāuatahanui and <1 % at Okatakata and Awanui.

Figure 7 Relative abundance plot of n-alkanes at Awanui, Okatakata and Pāuatahanui.

Total concentrations of C₂₁-C₃₃ n-alkanes in sediments ranged from 0.01 to 2.9 µg g⁻¹ TOC at Pāuatahanui, 0.01 to 4.4 µg g⁻¹ TOC at Okatakata, and from 0.7 to 9.4 µg g⁻¹ TOC at Awanui. Across sites, CPI ranged from 1.2 at Pāuatahanui to 10.8 at Awanui. OEP index was >2.2, ranging from 2.2 at Pāuatahanui to 51.4 at Awanui. ACL values varied between 27.2 and 29.9. P_aq ranged from 0.2 at Pāuatahanui to 1.0 at Awanui. Stanol–sterol ratios were generally low, ranging from 0 to 0.18, with deeper samples typically exhibiting higher values. C₂₈ steroids, specifically campesterol and its degradation product campestanol, and C₂₉ steroids, specifically sitosterol and its degradation product sitostanol, are used to calculate the stanol-sterol ratios, as these were the most abundant steroids in the samples. Stigmasterol was present in some samples but was too low in abundance to include in the ratio calculations. Table 3 summarises site-specific ranges, and Fig. S6 (Supplement) shows boxplots of biomarker indices for each core examined in the study. Example chromatograms showing the distribution of n-alkanes and steroids are presented in Figs. S7 and S8 (Supplement).

Table 3 Biomarker index ranges at Pāuatahanui, Okatakata and Awanui. Mean ± SE are presented in brackets.

4.8 Ramped-Pyrolysis Oxidation-Accelerator Mass Spectrometry (RPO-AMS) and Pyrolysis-Gas Chromatography-Mass Spectrometry (Py-GC-MS)

CRAs for the collected and measured pyrolytic splits ranged from 515±85 years Before Present (BP) to 2491±165 years BP at Okatakata, with measured pyrolysis temperatures ranging from 105 to 700 °C. At Awanui, CRAs ranged from 2350±160 years BP to 3246±188 years BP, with measured pyrolysis temperatures ranging from 105 to 700 °C. Isotopic mixing model results provide estimates of 67 % (first collected pyrolytic split) to 0 % (fifth collected pyrolytic split) and 62 % to 0 % of syndepositional OC within sediments at Okatakata and Awanui, respectively. Table 4 shows the RPO-AMS and isotopic mixing results, and Fig. S9 (Supplement) shows the RPO thermographs with the pyrolytic splits.

Figure 8 shows the identified relative abundances of the determined compound groups in all five measured pyrolytic splits. Most common compound types in the first split are polysaccharide derivatives, thiophenes and N-compounds, making up 72 % and 85 % at Okatakata and Awanui, respectively. Higher temperature splits show increasing proportions of recalcitrant OC (phenols, cyclic alkanes and alkylbenzenes) and undiagnostic sources. By the third split, these contribute 59 % and 48 % at Okatakata and Awanui, respectively. Figure S10 (Supplement) shows example chromatograms for the temperature splits following ramped Py-GC-MS analysis.

Table 4 RPO-AMS and isotopic mixing model (C_modern) results for Okatakata (Oka) and Awanui (Awa) samples.

Figure 8 Bar graphs of OC source composition for each pyrolytic split (°C) for (a) Okatakata sample Oka-MM1 (28–30 cm) and (b) Awanui sample Awa-MM2 (90–95 cm).

4.9 Principal Component Analysis (PCA) and hierarchical clustering

PCA of elemental, isotope, and XRF datasets for all samples (excluding δ¹⁵N, Mn, Sr, K and Ca due to their weak correlations) explained 61.3 % of the variance along PC1 and PC2 (Fig. 9a). PC1 reflects OM content with negative correlations for TOC, TN, CD, and positive correlations for δ¹³C_org, DBD, Si and Zr:Rb. TOC and TN displayed the highest loadings and Cos2 values, suggesting they are the most influential variables in PC1. PC2 reflects lithogenic content with positive correlations for Ti, Fe, Al, S, and δ¹³C_org. Ti and Fe had the highest loadings and Cos2 values, suggesting that they are the most influential variables in PC2. Hierarchical cluster analysis identified three distinct groups (Fig. 9b; Tables S9–S10). Cluster 1 (top 20 cm samples) exhibited statistically significantly (p<0.05) higher mean TOC, TN, and C:N, and lower δ¹³C_org, DBD, and Zr:Rb. Cluster 2 (0–50 cm) showed lower TOC, TN and C:N and higher DBD. Cluster 3 (5–95 cm) showed higher mean Al, Fe, Ti, Si, S, and δ¹³C_org.

PCA of elemental, isotope, XRF, and lipid biomarker datasets (excluding δ¹⁵N, Mn, Ca, Al and stanol-sterol due to their weak correlations) explained 64.3 % of the variance along PC1 and PC2 (Fig. 9c). PC1 represents lithogenic content with positive loadings for Ti, Fe, S, CPI, OEP, and P_aq, and negative correlations for Sr, Si, K, and ACL. PC2 is associated with OM content, showing negative correlations for TOC, TN, CD, C:N, and positive correlations for Zr:Rb, DBD and δ¹³C_org. Hierarchical clustering revealed three distinct groups with statistically significant (p<0.05) differences in means (Fig. 9d; Tables S11–S12). Cluster 1 (top 15 cm from Pāuatahanui, Okatakata, and Awanui) showed higher mean TOC and TN but lower DBD, Si, and Zr:Rb. Cluster 2 (0–45 cm from Pāuatahanui and Okatakata) had lower mean TOC, TN, S, Fe, CPI, and OEP, but higher Sr, Si, K and Zr:Rb, while Cluster 3 (5–95 cm Awanui and Okatakata) exhibited higher δ¹³C_org, S, Fe, Ti, OEP and P_aq, and lower ACL.

Figure 9 (a) PCA plots and (b) hierarchical cluster analysis of elemental, isotopic and XRF data sets for all study sites. (c) PCA plots and (d) hierarchical cluster analysis of elemental, isotopic, XRF and biomarker data sets. Clusters identified in hierarchical cluster analysis are outlined on the PCA plots.

5 Discussion

Data collected in this study provides important new information regarding carbon stocks and accumulation rates, and OM source and preservation at five saltmarsh locations in NZ. Soil OC stocks and CARs are highly variable, ranging from 41.3±9.4 to 92.3±66.2 Mg C ha⁻¹ and 0.46±0.02 to 1.53±0.09 Mg C ha⁻¹ yr⁻¹ (mean ± SE), respectively. These values are similar to previously reported values for saltmarshes and mangroves in NZ (Bulmer et al., 2024), but are smaller compared to global saltmarsh averages (Chmura et al., 2003; Maxwell et al., 2024; Ouyang and Lee, 2014; Wang et al., 2021). Vegetation type (zonation) appears to have no influence on TOC stock (Tables S2–S7), but accumulation of OC in the top 15–20 cm of sedimentary packages, proximity to allochthonous inputs, and site-specific geochemistry all appear to influence site variability (Fig. 9). Results from carbon and nitrogen stable isotope (Fig. 6) and lipid biomarker analyses (Fig. 7) indicate substantial contributions from saltmarsh vegetation to the OC pool. OM decay increases with depth at all sites (Fig. 9). However, abundant labile plant-derived OC in the lowermost saltmarsh sediments suggests that OC buried under anoxic conditions can be preserved over longer timeframes (Fig. 8). Figure 10 provides an overview of the key findings.

Figure 10 Schematic summary of key processes controlling organic carbon storage in saltmarsh sediments. The figure was created by the authors and synthesises concepts described by Howard et al. (2023).

5.1 Carbon stocks and accumulation rates

TOC stocks at all sites except Awanui are lower than estimated global representative soil OC stocks (83.1 Mg C ha⁻¹) derived from analysis of soils with an average depth of 30 cm (Maxwell et al., 2024). Low values at Robert Findlay are likely because saltmarsh deposits at this location are only 9 cm thick on average. Saltmarsh sediments at Okatakata and Pāuatahanui are reasonably thick (25 and 23 cm on average), and lower TOC stocks are likely due to site-specific geomorphic settings, such as lower siliciclastic inputs at Okatakata and historic anthropogenic impacts at Pāuatahanui. Saltmarsh deposits at Omaia are 23 cm thick on average, but were subject to increased OM decay rates due to greater exposure to oxygen when the marsh was converted to pasture (Ewers Lewis et al., 2018; Heuscher et al., 2005). Therefore, these TOC stocks likely reflect a higher-than-average amount of degraded marsh materials and soils that accumulated following drainage. Awanui recorded the highest mean TOC stock and the largest range in values, due to cores collected from juvenile to well-established marshes. Awanui also contains the thickest recovered sequence of saltmarsh sediments (95 cm). This deposit has a TOC stock value (421 Mg C ha⁻¹) that is almost five times higher than the national average for saltmarsh soils of similar thickness (92.5±12.42 Mg C ha⁻¹; Bulmer et al., 2024), and nearly double the estimated global representative stock for 1 m saltmarsh deposits (268 Mg C ha⁻¹; Maxwell et al., 2024).

Okatakata and Awanui saltmarshes are less than 700 m apart, but the mean TOC stock at Awanui is more than double that at Okatakata. Basal sediments at both sites have similar ages (1750±32 CE at Awanui and 1759±54 CE at Okatakata), but the mean sediment accumulation rate is 1.1 mm yr⁻¹ at Okatakata and 3.5 mm yr⁻¹ at Awanui. Awanui saltmarsh is closer to Awanui River than Okatakata and likely receives a relatively high input of allochthonous sediment. Previous studies have also attributed higher TOC stocks in saltmarshes that are adjacent to freshwater sources to the addition of terrestrial sediments (Hansen et al., 2017; Hayes et al., 2017; Kelleway et al., 2016; Peck et al., 2025; Van de Broek et al., 2016, 2018). Higher rates of mineral sediment deposition lead to more rapid burial of OC, which enhances its preservation (Kirwan and Mudd, 2012; Van de Broek et al., 2016, 2018) and contributes to the higher stocks (Fig. 10).

Mean CARs at our study sites vary between 0.46±0.02 to 1.53±0.09 Mg C ha⁻¹ yr⁻¹. These rates are lower than global average values but are comparable to the NZ mean estimate of 0.89 Mg C ha yr⁻¹ (Bulmer et al., 2024). A preliminary assessment of the restoration effect at Pāuatahanui and Robert Findlay shows increased carbon accumulation, with post-restoration rates approaching the highest observed values in each core (Figs. 10 and S11). Similar trends have been reported in restored saltmarshes in the United Kingdom and northwest Europe, where CARs increased by a factor of 1.6 following restoration (Drexler et al., 2020; Mason et al., 2022; Miller et al., 2022; Mossman et al., 2022).

Several factors can contribute to low CARs observed in saltmarshes around the world. For example, young saltmarsh deposits at locations in Great Britain (Smeaton et al., 2024) and Scandinavia (Leiva-Dueñas et al., 2024) have low mean (1.1±0.43 Mg C ha yr⁻¹) and median (0.32 Mg C ha yr⁻¹) accumulation rates, respectively. These low rates are attributed to the relatively young age of these marshes and the fact that relatively thin saltmarsh deposits form under temperate climatic conditions. Saltmarsh deposits in the southwest Atlantic (Martinetto et al., 2023) also have low mean carbon burial rates (0.48 Mg C ha⁻¹ yr⁻¹) that are likely caused by biological activity, including bioturbation by burrowing crabs that can mix and oxygenate the soils, causing OC degradation. The rate of relative sea level rise can also affect accumulation rates. Coastal sediments require accommodation space in which to accumulate, and this accommodation space is created as relative sea level rises (Rogers et al., 2019, 2022). Saltmarsh deposits in NZ seldom exceed 0.5 meters in thickness (Gehrels et al., 2008), partly because sea level has remained relatively stable across the Southern Hemisphere over the past 6 kyr (Rogers et al., 2023). All saltmarsh deposits examined in this study are young, bioturbated, and accumulated at locations where relative sea level has remained relatively stable. These factors all contribute to relatively low carbon stocks and accumulation rates.

5.2 Sources of organic matter in saltmarshes

Knowing the source of OM in saltmarsh soils can help us understand processes that drive carbon accumulation and can better inform the carbon mitigation potential of these ecosystems. δ¹³C_org and C:N data indicate that the OC accumulated in saltmarsh soils examined in this study is primarily derived from a combination of saltmarsh C₃ plants and estuarine biota (Fig. 6). Sediments at Okatakata, Omaia and Awanui also include OM from freshwater algae that was most likely transported to the sites via plumes from the Awanui River. Omaia sediments contain carbon derived from the widest range of sources, which likely reflects the influence of human activity, including drainage and farming, which mixed older marine sediments with saltmarsh deposits, younger topsoil and pasture vegetation. Although δ¹³C signatures presented in this study typically indicate C₃ plant inputs, these values cannot distinguish between saltmarsh vegetation, which is also C₃ in NZ, and other terrestrial C₃ plants, which include all indigenous NZ plants and at least some of the commercial exotic crops, such as rye grass Lolium perenne. Because all cores, except at Omaia, were collected within sites dominated by saltmarsh species (with surrounding mangroves at all sites except Pāuatahanui) and little/no adjacent forestry or terrestrial vegetation, we interpret the C₃ signal as primarily saltmarsh derived. However, we acknowledge that contributions from transported terrestrial material cannot be excluded, as tidal mixing and river discharge can deliver allochthonous marine- and terrigenous-derived OM (Alongi, 1997; Bulmer et al., 2020; Pondell and Canuel, 2022; Smeaton et al., 2024).

Given the overlap in observed isotopic and elemental ratios, it is challenging to quantify the relative contribution of allochthonous and autochthonous sources without knowing the specific δ¹³C_org, δ¹⁵N and C:N values for each end-member and their decomposition rates (Kumar et al., 2020b; Wang et al., 2003). These data were not collected in the study. However, lipid biomarker analysis offers additional insights into the origins of OM (e.g., Naeher et al., 2022; Peters et al., 2007). Odd-chain-length n-alkanes in the range C₂₁ to C₃₃ dominate in Pāuatahanui (98 %), Okatakata (99 %), and Awanui (99 %) and indicate epicuticular waxes from saltmarsh plants (Fig. 7). Biomarkers from marine and freshwater algae, bacteria and/or other microorganisms (< C₂₀) (Cranwell, 1981; Eglinton and Hamilton, 1967; Ficken et al., 2000; Kumar et al., 2019; Meyers, 1997; Zhang et al., 2024) were insignificant across all sites. Positive correlations between CPI and OEP with values >3 in surface sediments also suggest that saltmarsh vegetation contributes to the sedimentary OM (Figs. 9 and 10; Kennicutt et al., 1987; Zhao et al., 2024).

Biomarkers can also help identify which individual saltmarsh plant species the OM originates from. For example, the dominance of C₂₃ and C₂₅ n-alkanes and P_aq values >0.6 at Okatakata and Awanui suggest soil organic material is derived from Salicornia quinqueflora (Tanner et al., 2007, 2010). Other aquatic plants that could be contributing to the mid-chain n-alkanes include moss species Kindbergia spp. and Polytrichum spp. (Ortiz et al., 2016), both of which were observed in low and mid marsh areas. It is also possible that floating macrophytes are delivered via tidal/riverine inputs; however, these were not observed during data collection. In contrast, dominant peaks at C₂₇, C₂₉, C₃₁ and C₃₃, and P_aq<0.6 at Pāuatahanui reflect emergent macrophytes and salt-tolerant terrestrial plants. Studies have demonstrated that C₃ saltmarsh grass, rush, sedge and shrub species exhibit maximum abundances in the C₂₇–C₃₃ range, consistent with the n-alkane distributions, as well as ACL values of 27–29 (Eley et al., 2016; Ferreira et al., 2009; Ortiz et al., 2011; Pondell and Canuel, 2022; Tanner et al., 2007, 2010; Wang et al., 2003; Zhang and Wang, 2019). PCA plots confirm these observations as ACL and P_aq display opposite trends, and an increase in P_aq and a decrease in ACL indicate dominance of non-emergent macrophytes (Fig. 9). Upstream terrestrial plants, such as Pinus radiata can also be a source of C₂₉ n-alkane (Gonzalez-Vila et al., 2003; Kumar et al., 2020a). However, there are no pine plantations adjacent to the study sites, and terrestrial plant debris/litter was not observed during fieldwork. Seagrass species may contribute C₂₉ n-alkanes to the OC pool (Jaffé et al., 2001; Kumar et al., 2020a). Seagrass habitats are present at Pāuatahanui (Zabarte-Maeztu et al., 2020), and Rangaunu Harbour has one of the largest seagrass habitats in the country (Bulmer et al., 2024). Therefore, cross-habitat exchange could contribute to the composition of the OC pool (Bulmer et al., 2020).

Furthermore, plant-derived OC tends to be more resistant to microbial breakdown than OC from algae and bacteria, which may result in its preferential preservation in sediments (Schmidt et al., 2011; Zhang et al., 2024). Clearly, the estimated CARs presented in this study reflect both autochthonous and allochthonous carbon input and do not reflect the amount of OC directly sequestered from the atmosphere through in-situ production (Smeaton et al., 2024). However, lipid biomarker analysis demonstrates that saltmarsh vegetation makes up a substantial portion of OM at the studied sites, and burial and preservation of autochthonous OM contribute to the long-term carbon storage capacity of saltmarshes.

5.3 Organic matter preservation

Strong positive correlation between TOC and TN at all sites suggests they are likely influenced by similar geochemical processes (Fig. 9; Brandini et al., 2022; Zhang et al., 2024). Surface samples generally exhibit higher mean OM contents (TOC, TN, CD, C:N) and lower mean concentrations of lithogenic components (Si, K, Ti, Fe, Zr:Rb). This indicates that young, finer-grained, organic-rich surface sediments with fresh OC are predominant across study sites (Kelleway et al., 2017; Krüger et al., 2024; Mazarrasa et al., 2023). Finer-grained samples typically have higher OC content because a higher proportion of silt and clay in sediments enhances preservation of OM (Mazarrasa et al., 2023; Meyers, 1994; Russell et al., 2023). In contrast, older, deeper, coarser-grained samples generally show lower OM, higher lithogenic content, enriched δ¹³C_org, higher stanol-sterol values, lower C:N ratios, and lower CPI values of 1–2. These characteristics suggest that OM has been preferentially utilised and decomposed by microbial activity over time, increasing the lithogenic contribution in the OC pool (Benner et al., 1987; Jaffé et al., 2001; Krüger et al., 2024; Zhao et al., 2024). Omaia shows enriched δ¹³C_org, low TOC, TN, and CD, and increased DBD values in the majority of samples, reflecting the degraded and compacted nature of the soils.

At Awanui and Okatakata, deeper samples display higher OM content alongside sulphur, suggesting preservation of OC under sulphur-rich anoxic conditions (Antler et al., 2019; Froelich et al., 1979; Thamdrup et al., 1994). This aligns with the observation that oxic and suboxic degradation of OC is restricted to shallow depths (<20 cm), and deeper soil horizons remain more consistently anaerobic, promoting preservation (Howarth and Teal, 1979; Spivak et al., 2019). However, the dominance of aerobic or anaerobic conditions also depends on the position of the saltmarsh in the tidal frame. Locations that are flooded less frequently tend to oxidise to depths up to 30 cm below the surface (e.g., Mueller et al., 2019). In these “drier” locations, other factors, such as microbial community compositions, formation of organo-mineral complexes, and physical protection by mineral aggregates, may control OM preservation (Barber et al., 2017; Macreadie et al., 2025; Spivak et al., 2019). Samples recovered from lower depths in the cores taken at Awanui and Okatakata also have the highest Ti, Al, and Fe contents. Our age-depth models indicate these sediments were deposited when flood control and drainage works were carried out in the Awanui River catchment between 1916 and 1936. These activities likely increased downstream transport of terrestrial sediment to the harbour (Cathcart, 2005). It is possible that deposition of mineral sediments and their geochemical interactions with OC (e.g., OC binding with iron oxides to form stable Fe-OC complexes) enhanced OM preservation (Barber et al., 2017; Macreadie et al., 2025).

Isotopic mixing model results from RPO-AMS analysis provide a first-order approximation of the relative contributions of labile, syndepositional versus older, recalcitrant OC in the soils. The results suggest a relative proportion of two-thirds labile and one-third recalcitrant OC at Okatakata and Awanui basal samples for the lowest-temperature split (Table 4). RPO-AMS split age profiles of both basal samples at Okatakata and Awanui show a wide range of ages that are older than those estimated using ²¹⁰Pb age-depth models, representing the mixture of labile and recalcitrant carbon delivered to the site. Awanui has much older carbon compared to Okatakata. This is consistent with the core's proximity to the river, which delivers allochthonous material from the upper catchment.

Py-GC-MS results provide additional insights into the OM composition in samples that were also processed for RPO-AMS ages. Lower-temperature splits generated from basal samples from Okatakata (100 to 300 and 300 to 385 °C) and Awanui (100 to 280 and 280 to 333 °C) show higher relative contributions of syndepositional OC. These lower-temperature splits release more labile and volatile compounds, which are more closely associated with the sources and state of syndepositional carbon (Ginnane et al., 2024; Rosenheim et al., 2008). For example, furans indicate the presence of higher-molecular-weight polysaccharides such as cellulose, a major carbohydrate constituent of plant cell walls that is readily consumed by microbes (Carr et al., 2010; Kaal et al., 2020). Further heating during pyrolysis (typically >400 °C) breaks down the macromolecular structure of OM compounds, releasing more recalcitrant and older reworked OM (Ginnane et al., 2024; Maier et al., 2025). Phenols, for instance, reflect lignin, a macromolecular compound found almost exclusively in the cell walls of terrestrial vascular plants, which is more resistant to OM degradation (Grandy and Neff, 2008; Kaal et al., 2020; Zhang and Wang, 2019). Higher-temperature pyrolytic splits for Okatakata and Awanui also produce higher relative abundances of cyclic alkanes/alkylbenzenes and other undiagnostic or refractory compounds (Fig. 8). Collectively, Py-GC-MS results indicate preservation of both labile and recalcitrant saltmarsh plant-derived OC under anoxic conditions at Okatakata and Awanui (Figs. 8 and 10; González-Pérez et al., 2012). However, N-containing pyrolysis products (e.g., indole, benzonitrile, pyridine) in coastal wetland soils have also been attributed to amino acids and proteins from algae and phytoplankton (Carr et al., 2010; Kaal et al., 2020), as well as bacterial biomass within soil OM (Ferreira et al., 2009; Zhang and Wang, 2019; Zhu et al., 2016). TLE and Py-GC-MS analyses of the basal samples suggest ≤1 % and ≤5 % contributions from marine source n-alkanes at Okatakata and Awanui, respectively (Figs. 7 and 8), indicating that N-compounds in the samples likely correspond to microbial biomass. These findings suggest some degree of post-depositional microbial decomposition, but in-situ plant biomass preserved under anoxic conditions is evident in the signal.

5.4 Implications and future research

This study provides the first integrated assessment of carbon stocks, accumulation rates, and OM source and preservation across a latitudinal gradient of saltmarshes in NZ. Our findings reveal that carbon accumulation and preservation are strongly influenced by site-specific factors such as land use history, sediment characteristics, and allochthonous inputs. Notably, we demonstrate that even relatively young saltmarsh deposits can store substantial amounts of carbon when conditions favour OM preservation. Our results also corroborate previous research showing that soil OM properties, carbon stocks and accumulation rates vary widely with geomorphic setting, land use history, tidal regime, salinity, vegetation type, local lithologic input, and soil depth (e.g., Ewers Lewis et al., 2019; Hansen et al., 2017; Kelleway et al., 2016; Martinetto et al., 2023; Martins et al., 2022; McMahon et al., 2023; Owers et al., 2020; Ruiz-Fernández et al., 2018; Russell et al., 2023; Saintilan et al., 2013). Blue carbon sampling strategies must consider these sources of variability to achieve accurate assessment. Our results suggest that allochthonous inputs, particularly in locations proximal to freshwater sources such as rivers, play a critical role in enhancing carbon accumulation and preservation. The study also highlights the importance of measuring OC to the base of the saltmarsh deposit to capture the full extent of carbon storage and preservation processes. Future research could explore the role of other environmental variables such as distance to tidal creeks and freshwater sources, vegetation biomass, elevation, inundation frequency and duration, salinity and sediment composition on OM dynamics and carbon stocks in NZ saltmarshes (e.g., Hansen et al., 2017; Janousek et al., 2025; McMahon et al., 2023; Puppin et al., 2024; Russell et al., 2023).

These findings have practical implications for blue carbon assessment methodologies. Current blue carbon methodologies, such as those used by Verra, often require removal of allochthonous material to estimate accumulation rates (Needelman et al., 2018). Results from this and other recent studies (e.g., Li et al., 2025; Peck et al., 2025) suggest that allochthonous inputs can enhance long-term preservation of OC. Furthermore, allochthonous carbon transported to restored wetlands in tidally influenced locations is considered “additional” because its preservation results directly from the restoration activities (Lovelock et al., 2023a, b). Strict exclusion of allochthonous contribution to saltmarsh OM, when these inputs are not already accounted elsewhere, may underestimate the true sequestration potential of saltmarsh ecosystems. Streamlined protocols that balance scientific rigour with practical requirements are needed to improve carbon mitigation estimates and scale BCE restoration globally.

6 Conclusions

This study successfully quantified carbon stocks and accumulation rates across five saltmarsh sites in NZ, revealing significant variability influenced by site-specific factors such as geomorphic settings, land use, and proximity to allochthonous inputs. Mean carbon stocks ranged from 41.3±9.4 to 92.3±66.2 Mg C ha⁻¹, and mean accumulation rates ranged from 0.46±0.02 to 1.53±0.09 Mg C ha⁻¹ yr⁻¹ (mean ± SE). Although these values are lower than global averages, they are comparable to previous national estimates for NZ. Stable isotopes, lipid biomarkers, and RPO-AMS combined with Py-GC-MS analyses provided complementary insights into the source and composition of soil OM. Together, these analyses indicate that in-situ saltmarsh vegetation is a major contributor to the OM pool and that plant-derived OC remains preserved in saltmarsh soils for several centuries. Overall, these findings improve national estimates of carbon accumulation in saltmarsh ecosystems and advance the methodological approaches needed to assess the potential for BCEs to capture carbon and mitigate climate change impacts.