Quantification of blue carbon in salt marshes of the Pacific coast of Canada

Tidal salt marshes are known to accumulate “blue carbon” at high rates relative to their surface area, which render these systems among the Earth's most efficient carbon (C) sinks. However, the potential for tidal salt marshes to mitigate global warming remains poorly constrained because of the lack of representative sampling of tidal marshes from around the globe, inadequate areal extent estimations, and inappropriate dating methods for accurately estimating C accumulation rates. Here we provide the first estimates of organic C storage and accumulation rates in salt marshes along the Pacific coast of Canada, within the United Nations Educational, Scientific and Cultural Organization (UNESCO) Clayoquot Sound Biosphere Reserve and Pacific Rim National Park Reserve, a region currently underrepresented in global compilations. Within the context of other sites from the Pacific coast of North America, these young Clayoquot Sound marshes have relatively low C stocks but are accumulating C at rates that are higher than the global average with pronounced differences between high and low marsh habitats. The average C stock calculated during the past 30 years is 54 ± 5 Mg C ha⁻¹ (mean ± standard error), which accounts for 81 % of the C accumulated to the base of the marsh peat layer (67 ± 9 Mg C ha⁻¹). The total C stock is just under one-third of previous global estimates of salt marsh C stocks, likely due to the shallow depth and young age of the marsh. In contrast, the average C accumulation rate (CAR) (184 ± 50 g C m⁻² yr⁻¹ to the base of the peat layer) is higher than both CARs from salt marshes along the Pacific coast (112 ± 12 g C m⁻² yr⁻¹) and global estimates (91 ± 7 g C m⁻² yr⁻¹). This difference was even more pronounced when we considered individual marsh zones: CARs were significantly greater in high marsh (303 ± 45 g C m⁻² yr⁻¹) compared to the low marsh sediments (63 ± 6 g C m⁻² yr⁻¹), an observation unique to Clayoquot Sound among NE Pacific coast marsh studies. We attribute low CARs in the low marsh zones to shallow-rooting vegetation, reduced terrestrial sediment inputs, negative relative sea level rise in the region, and enhanced erosional processes. Per hectare, CARs in Clayoquot Sound marsh soils are approximately 2–7 times greater than C uptake rates based on net ecosystem productivity in Canadian boreal forests, which highlights their potential importance as C reservoirs and the need to consider their C accumulation capacity as a climate mitigation co-benefit when conserving for other salt marsh ecosystem services.

1 Introduction

Coastal, vegetated ecosystems, such as seagrass meadows, mangroves, and tidal salt marshes, have recently been recognized for their ability to store large amounts of “blue carbon” within their soils and sediments (Kennedy et al., 2013; Howard et al., 2014). While blue carbon ecosystems cover approximately 0.2 % of the ocean surface, they have been estimated to be responsible for up to 50 % of total coastal ocean carbon (C) burial (including estuaries and continental platforms) (Duarte et al., 2005), and their per-area C sequestration rate is substantially greater than that of terrestrial forest soils (McLeod et al., 2011). Globally, blue C ecosystems have been estimated to sequester between 75 and 224 Tg C yr⁻¹ (Duarte et al., 2013). For comparison, deep ocean organic carbon burial rates were recently estimated to be 20 ± 6 Tg C yr⁻¹ (Hayes et al., 2021), with an additional burial of 222 Tg C yr⁻¹ estimated for shelves and platforms shallower than 1000 m (Burdige, 2007). Due to this high C storage and accumulation rate capacity per unit area, coastal vegetated ecosystems have been suggested to play an important role in climate warming mitigation (Howard et al., 2017). However, when the ecosystem is degraded, the stored C can be released, and annual C uptake by the ecosystem ceases, resulting in losses of ecosystem services (McLeod et al., 2011; Pendleton et al., 2012). Thus, to better inform policies that identify priority areas for conservation, more precise measurement of C stocks and accumulation potential are needed (Howard et al., 2017).

Global estimates of salt marsh area, C stocks, and C accumulation rates (CAR) are subject to large uncertainties. Duarte et al. (2013) noted a 20-fold uncertainty in global estimates of salt marsh area (ranging from 22 000 to 400 000 km²) associated with ambiguous classification schemes for wetlands. For example, some classification systems consider freshwater and saltwater marshes in the same category (Duarte et al., 2013). Similarly, the estimated global soil C stock of all salt marshes ranges between 0.4 and 6.5 Pg C, a 16-fold range (Duarte et al., 2013). The average global CAR for salt marshes was estimated as 91 ± 19 g C m⁻² yr⁻¹ by the Intergovernmental Panel on Climate Change (IPCC) (Kennedy et al., 2013) and then as 245 ± 26 g C m⁻² yr⁻¹ in a subsequent global compilation (mean ± standard deviation) (Ouyang and Lee, 2014). Furthermore, these reviews of salt marsh CAR estimates note disproportionate representation from certain areas of the world (Ouyang and Lee, 2014; Chmura et al., 2003). Some areas, such as Europe and eastern North America, have dozens of CAR data points, while others, such as western North America, East Asia, and Australia are underrepresented. Regions such as Africa, India, and South America have fewer or no data at all. The high variability in CAR from site to site combined with the 20-fold uncertainty in global marsh area estimates result in global salt marsh CAR estimates ranging from 0.9 to 31.4 Tg C yr⁻¹ (Ouyang and Lee, 2014). This 35-fold range is 7 times greater than the global range for mangroves (Ouyang and Lee, 2014; Donato et al., 2011).

The extensive use of ¹³⁷Cs radioisotope or a marker horizon method for sediment dating also limits comparisons of CAR between studies. The difficulties of quantifying low concentrations of ¹³⁷Cs in coastal sediments, and its low retention in high-organic-content sediments (Davis et al., 1984), limit the applicability of this technique to estimate CAR in salt marshes and lead to overestimates of sediment accumulation rates when compared to other dating methods such as ²¹⁰Pb dating (e.g. Callaway et al., 2012; Johannessen and MacDonald, 2016). Unlike ¹³⁷Cs, the natural radionuclide ²¹⁰Pb is highly retained in organic matter, which makes it suitable for establishing the chronology of sedimentary deposits accumulated over approximately the past 100 years (Krishnaswamy et al., 1971; Arias-Ortiz et al., 2018). Also, the lower detection limits of ²¹⁰Pb reduce analytical errors in the derived chronologies below those obtained by using ¹³⁷Cs thus providing better constrained recent CAR for salt marshes (e.g. Corbett and Walsh, 2015).

The Commission for Environmental Cooperation (CEC), a tri-national governmental organization promoting scientific cooperation between Canada, the United States, and Mexico, identified the Pacific coast of Canada as a significant data gap in terms of our knowledge of the aerial extent and quantification of blue C storage and accumulation. Additionally, a review of global salt marsh CAR data identified only eight sites on the entire Pacific coast of the continent, none of which were north of 38.2^∘ N (Ouyang and Lee, 2014). This study aims to address this data gap by providing C stock and CAR from the Pacific coast of Canada as a part of the government of Canada's contribution to a continent-wide assessment of blue C mitigation potential. We sampled seven salt marshes within the United Nations Educational, Scientific, and Cultural Organization (UNESCO) Clayoquot Sound Biosphere Reserve; British Columbia's Tofino Mudflats Wildlife Management Area; and Pacific Rim National Park Reserve of Canada on Vancouver Island, British Columbia (Fig. 1). These mesotidal estuarine marshes, often constrained in size by surrounding topography, are typical of the marshes found on the Pacific coast of British Columbia and therefore provide a good representation of many of the tidal wetland ecosystems found here (Ryder et al., 2007). We calculated soil C stocks (per unit area) across the high and low marsh zones of each marsh and used ²¹⁰Pb dating in a subset of the collected cores to quantify CAR. We then used the extent of high and low marsh areas from aerial imagery to estimate total C storage and total annual C accumulation. Finally, we placed these new data within the context of CAR estimates from salt marshes in the NE Pacific region.

Figure 1 Site locations. Cypress River Flats (CRF), Cannery Bay East (CBE), Cannery Bay West (CBW), Kennedy Cove South (KCS), Shipwreck Cove (SWC), Grice Bay at Kootowis Creek (GBK), and Tofino Mudflats (TMF). Pacific Rim National Park and Reserve (blue crosshatching) covers the southern portion of the map and the Tofino Mudflats Wildlife Management Area (pink crosshatching) covers portions of the southwestern area. The study region lies within the UNESCO Clayoquot Sound Biosphere Reserve, British Columbia, Canada (inset, outlined).

2 Methods

2.1 Study area

Clayoquot Sound is a complex system of inlets on the west coast of Vancouver Island, British Columbia, Canada. The sparsely populated area includes several protected area designations, including the Long Beach Unit of Pacific Rim National Park and Reserve of Canada, the Province of British Columbia's Tofino Mudflats Wildlife Management Area, and the UNESCO Clayoquot Sound Biosphere Reserve, which protects 366 000 ha of the west coast of Vancouver Island (Fig. 1). The region is part of the temperate rainforest biome with high annual rainfall (3270 mm yr⁻¹) and average annual temperature of 9.5 ^∘C (Environment Canada, 1981–2010 averages for Tofino). While we were not able to measure the precise tidal range at each location, the mean tidal range in Tofino is approximately 2.7 m (DFO 2022).

These sites are typical of salt marshes along Canada's Pacific coast because they include small marshes along protected shorelines and bays as well as larger estuarine marshes near creeks and rivers. All sites were somewhat close to fluvial sources of varying size. Surface water salinity in the surrounding waters ranged from 5.9 at KCS to 24 in Grice Bay and 29 at Roberts Point 6 km south of CRF (Postlethwaite et al., 2018).

2.2 Field sampling

Within each marsh (n=7), vegetation composition was noted in 50 × 50 cm quadrats before sediment cores (n=34) were extracted during summer 2016 along linear transects perpendicular to the low-tide shoreline, following the methodology of Howard et al. (2014). Coring spots were approximately evenly spaced along the transect (between 9 and 24 m apart) from land to sea and spanned the low and high marsh zones. Core locations were chosen to avoid ditches and channels without organic soil accumulation; these ditches and channels made up a relatively small portion of total marsh surface area (<5 %).

We note that, within this system, the partitioning between high and low marsh zones appears very closely linked with associated vegetation. Therefore, vegetation composition was recorded as an indicator of “low” vs. “high” marsh zones (Porter, 1982; Weinmann et al., 1984) around each coring spot. Coring spots were considered low marsh if the species Triglochin maritima, Salicornia spp., Fucus ssp., or Distichlis spicata were present and high marsh if it included Plantago maritima, Deschampsia caespitosa, Grindelia integrifolia, Potentilla anserina, Lysimachia maritima, or Eleocharis ssp. If a spot contained a mixture of these species, the majority percent cover of high or low marsh species was used to determine whether the spot was high or low marsh. Carex lyngbyei were often found throughout both strata and so were not considered unique to one zone. The high marsh species' ranges align approximately with the mean extreme high-water line of estuarine marshes in Clayoquot Sound, while low marsh encompasses elevations between the mean lower high water and the mean extreme high-water lines (Jefferson, 1973, as cited in Deur, 2000; Weinmann et al., 1984). This method was verified using detrended correspondence analysis, which showed that vegetation assemblages reflected distinct low and high marsh zones (Hill and Gauch, 1980; Appendix B). This tight coupling between vegetation type and marsh zone has also been observed in other studies on the west coast of Vancouver Island (Deur et al., 2000) and studies in nearby Boundary Bay (Gailis et al., 2021). Vegetation associations related to salinity and inundation are well documented and commonly used in salt marsh delineation (MacKenzie and Moran, 2004).

Sediment cores were collected using a simple percussion coring technique in which a length of two in. (57 mm) diameter PVC vacuum tubing fitted with a plastic core catcher (AMS Inc.) was hammered into the ground until the depth of refusal (DoR). The DoR was reached between depths of 7 and 62 cm. In 19 of the 34 cores, the DoR was reached when the corer penetrated approximately 5–15 cm into layers of sand or gravel at the base of the marshes (see Appendix, Table A2). A similar layer is reached in the eelgrass meadows that form seaward of these marshes (Postlethwaite et al., 2018). At these locations, our cores were able to sample the depth of initiation of marsh-related organic C accumulation. At 8 of the 34 locations, the cores reached depths where peat mixed substantially with sand or clay; the remaining seven cores ended in peat layers. With three exceptions, the DoR occurred at minimum %C.

At one site (GBK) a steel sledge corer (AMS Inc.) was used to extract four cores, but mechanical problems required switching to the simpler method described above. These four cores penetrated approximately 15–30 cm into the sand/gravel layer below the marsh, which contained minimal amounts of C. Marsh characteristics (e.g. average percent carbon, dry bulk density, soil C density, C stocks, and C accumulation rates) were determined by considering only the upper part of the profile containing peat and overlying sand, gravel, or clay, and these cores were excluded from our average estimates of the DoR.

While in the field, all cores were stored upright in portable coolers with ice packs until their return to the laboratory where they were photographed, logged, and stored long-term under refrigeration at 4 ^∘C at the Parks Canada laboratory in Vancouver, British Columbia. Carbon and ²¹⁰Pb laboratory work was completed 2–4 months after sample collection, between July and October 2016. Additional dating laboratory work was completed by approximately 30 months after collection.

Use of the percussion corer resulted in sediment compaction during sample collection, which averaged about 20 % across all cores (range of 0 %–55 %) (Table A1 in the Appendix). Nevertheless, we opted to use a percussion corer instead of a gouge corer because the percussion corer had a closed chamber with internal PVC sleeves. Our experience with these sediments has shown that a gouge corer would have been susceptible to disturbance and sediment mixing due to the nature of the open chamber of the corer. We also did not use a Russian corer because compaction would have been similar to what we experienced with the percussion corer, and we did not want to introduce increased contamination through the pivoting nature of the sampling chamber with the Russian corer. Digging pits with a shovel was not an option as this study took place in a national park and biosphere reserve. We note below that correction for compaction was not necessary for estimation of C stocks because the C stocks were estimated directly from sediment cores and not from the overall depth of marsh soils (thus all carbon in the peat layer, regardless of compression, is included in the calculation). Furthermore, when we have estimated ²¹⁰Pb-derived accumulation rates (Fig. 6) we have done so in terms of cumulative mass (g cm⁻²) instead of depth (e.g. Gifford and Roderick, 2003). When we do need to account for compaction (e.g. Fig. 3), we use a compaction factor as described in Gailis et al. (2021), estimated for each core by dividing the length of core penetration by the length of core recovered (Table A1).

2.3 Marsh area estimations

ArcMap 10.3 tools were used with 50×50 cm resolution aerial orthophotos taken in July 2014 (Government of British Columbia) to estimate the area of high and low marsh zones. The difference between high marsh and low marsh was delineated by eye between darker-green, denser high marsh vegetation and lighter-green, salt-tolerant, and less-dense low marsh vegetation. This method was verified using the detrended correspondence analysis (e.g. Hill and Gauch, 1980) of vegetation survey data and was found to accurately categorize 94 % of the cores into the correct marsh zone (see Appendix B). Similar approaches in Boundary Bay marsh in BC, Canada, have resulted in a high degree of correspondence between marsh colour and low vs. high marsh delineation (Gailis et al., 2021).

2.4 Soil carbon content

Organic C content (%C = g C g⁻¹ marsh sediment) was determined using loss on ignition (%LOI) validated with carbon–nitrogen (CN) elemental and coulometric analysis (Froelich, 1980). LOI analysis was performed on every 1 cm subsample by homogenizing samples with a mortar and pestle, combusting them at 550 ^∘C for 4 h, weighing, and combusting again at 1000 ^∘C for 2 h (Heiri et al., 2001). The percentage mass loss-on-ignition (%LOI) was estimated as:

$$\begin{array}{}\text{(1)}& \mathit{\%}\text{LOI}={\displaystyle \frac{{\mathrm{DW}}_{\mathrm{i}}-{\mathrm{DW}}_{\mathrm{f}}}{{\mathrm{DW}}_{\mathrm{i}}}}\times \mathrm{100},\end{array}$$

where DW_i is initial dry weight (g) and DW_f is the dry weight (g) after burning. For comparison, the %C was also estimated by measuring total C (%TC) and inorganic C (%IC) on a subset of 93 samples (see Appendix A). %TC was measured on these homogenized subsamples using dry combustion elemental analysis with an Elementar elemental analyzer for CN analysis at the University of British Columbia's Department of Earth, Ocean and Atmospheric Sciences. The same subsamples were then analyzed for %IC using a UIC CM5014 CO₂ coulometer connected to a UIC CM5130 acidification module in the Climate, Oceans, and Paleo-Environments (COPE) laboratory at Simon Fraser University. Measurements of %IC were subtracted from the %TC measurements to estimate %C (Hodgson and Spooner, 2016; Hedges and Stern, 1984; Schumacher, 2002; Howard et al., 2014). Inorganic C was negligible in all 93 of the subsamples analyzed (max: 0.015 %) and assumed to be zero for all C calculation purposes. The strong correlation (R²=0.96, p<0.05) observed between %LOI and %C for this set of samples allowed us to use this relationship to convert %LOI to %C for all samples using the following equation, setting any negative %C value resulting from the use of a negative intercept equal to zero:

$$\begin{array}{}\text{(2)}& \mathit{\%}\text{C}=\mathrm{0.44}\left(\mathit{\%}\text{LOI}\right)-\mathrm{1.80}.\end{array}$$

Soil C density (SCD (g C cm⁻³)) was obtained as the product of the C content and the dry bulk density (DBD (g C cm⁻³)):

$$\begin{array}{}\text{(3)}& \text{SCD}\left({\displaystyle \frac{\mathrm{g}\mathrm{C}}{{\text{cm}}^{\mathrm{3}}}}\right)=\left({\displaystyle \frac{\mathit{\%}\text{C}}{\mathrm{100}}}\right)\times \text{DBD},\end{array}$$

where DBD was estimated after drying the sediment for no less than 72 h at 60 ^∘C. Specifically, DBD was measured by bisecting the core while still inside its casing with a clean, serrated knife to avoid cross-contamination or inaccurate volume measurement that might be caused by extruding the core. Volume was precisely measured immediately afterward using a custom-designed, brass implement designed to remove a single cubic centimeter of soil at a time from the flat face of the bisected core.

2.5 Marsh profiles chronology

The naturally occurring radionuclide ²¹⁰Pb was analyzed in sediment from a subset of eight cores to determine recent C accumulation rates (CARs) (Appleby et al., 1992). One high marsh and low marsh core each from CBE, CRF, GBK, and TMF were chosen to represent the two marsh zones in various sites of different sizes. We used these vegetation communities to establish representativeness of the high and low marsh sites that we sampled for ²¹⁰Pb analysis. A set of subsamples (n=6–7 for cores shorter than 20 cm, n=11–17 for longer cores) was sent to Core Scientific International (Winnipeg, Canada), Flett Research (Winnipeg, Canada), and MyCore Scientific (Dunrobin, Canada) to determine ²¹⁰Pb content. Activities of ²¹⁰Pb were determined by α spectrometry through its granddaughter ²¹⁰Po, assumed in secular equilibrium. The atmospheric or excess ²¹⁰Pb (²¹⁰Pb_xs) fraction used to derive the age–depth model was determined as the difference between the total ²¹⁰Pb activity and its parent nuclide ²²⁶Ra activity (Table A2). ²²⁶Ra activities were determined by α-spectrometry at Flett Research using calibrated geometries in a glass vessel, Spectech UCS 30 Alpha Scintillation Spectrometer purged with helium, and sealed for at least 11 d (Mathieu et al., 1988). Samples were sealed and stored for 2 h before counting to ensure secular equilibrium of ²²⁶Ra daughters. ²²⁶Ra was determined through counting ²²²Rn activity for 60 000 s (minimum detectable activity (MDA) = 0.0167 Bq kg⁻¹).

Because most of the cores did not show complete excess ²¹⁰Pb profile (see Results), sediment chronologies were estimated using the constant flux : constant sedimentation (CF : CS) model (Krishnaswamy et al., 1971; Eq. 4), which assumes a constant atmospheric deposition of ²¹⁰Pb to the marsh surface and a constant mass accumulation rate (MAR (g cm⁻² yr⁻¹)). Under these assumptions, the MAR can be obtained from the slope (cm² g⁻¹) of the linear best-fit line of the relationship between the natural log of the excess ²¹⁰Pb activity against the cumulative mass (g cm⁻²) using

$$\begin{array}{}\text{(4)}& \text{MAR}\left(\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\right)={\displaystyle \frac{\mathrm{0.0311}\left({\text{yr}}^{-\mathrm{1}}\right)}{\text{slope}\left({\text{cm}}^{\mathrm{2}}{\text{g}}^{-\mathrm{1}}\right)}},\end{array}$$

where 0.0311 yr⁻¹ is the radioactive decay constant of ²¹⁰Pb.

2.6 Carbon stocks

Organic C stocks (g cm⁻²) are used as a proxy for total organic C storage in marshes and are calculated by integrating all peat layers (n) of the profile (cm):

$$\begin{array}{}\text{(5)}& \text{C}\text{stock}\left(\mathrm{g}\mathrm{C}{\mathrm{cm}}^{-\mathrm{2}}\right)={\sum }_{i=\mathrm{0}}^n{\text{SCD}}_{\text{i}}\times \mathrm{1}\text{cm}.\end{array}$$

As most cores penetrated to a sand or gravel layer found at the base of the marsh (see Sect. 2.2), we are confident that the estimated C stocks (“base of peat C stocks”, hereafter) capture all marsh-associated C emplaced since the formation of the marsh. Because C stocks are estimated directly from sediment cores and not from the overall depth of marsh soils, correction of core compaction was not necessary.

To facilitate comparison with other studies, we also report the mean C stock (Mg C ha⁻¹) to a depth of 20 cm for the region (20 cm C stocks, hereafter). Here we estimated the accumulated C to the corrected (uncompacted) depth of 20 cm in each core, with the caveat that 3 of the 34 cores did not reach a depth of 20 cm and therefore were not included in this calculation (Table A1).

Measurements of C stocks based on the accumulation of organic C up to certain depth, however, might give an incorrect impression of C sequestration capacity when comparing sites with different accumulation rates. The dating of the cores, however, allows one to circumvent this problem by normalizing C stocks to a certain age horizon. We therefore used a horizon of 30 years before time of sampling to calculate these stocks (hereafter 30-year C stock). We chose this age–depth horizon because it represented the oldest ²¹⁰Pb date shared between the eight dated cores.

Finally, we estimated the total C storage (Mg C), in each marsh and for all seven marshes combined, by multiplying the average C stocks by the aerial extent of each marsh (Table 1). The total C storage was also estimated for low marsh and high marsh areas separately (Mg C). These estimates were made both to the depth of the base of the deepest peat (base of peat) layer and for the 20 cm horizon carbon stock estimates. Note that total C storage was not estimated for the 30-year horizon because we only had ages in four of the seven marshes, and thus there were three marshes where C storage could not be estimated to 30 years.

2.7 Carbon accumulation rates

We calculated C accumulation rates (CARs) in two ways: (1) to the “base of peat” depth of the deepest peat-containing layer overlying sand, gravel, or clay and (2) to the greatest, common age horizon shared by all cores to allow for comparisons between cores over equivalent time spans (30 years), which we call the “30-year CAR”. We include the whole-core CAR calculation to assess any potential bias in the results due to higher overall C content in more recent accumulation and for easier comparison with other studies. This organic C fraction is calculated as

$$\begin{array}{}\text{(6)}& \mathrm{\%}\mathrm{C}={\displaystyle \frac{\sum _{t=\mathrm{0}}^{t=i\text{years}}\text{C}\left(\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}\right)}{\sum _{t=\mathrm{0}}^{t=i\text{years}}\text{m}\left(\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}\right)}},\end{array}$$

where i is either 30 years or the oldest peat-containing layer with detectable excess ²¹⁰Pb activity for the core. The CAR is the product of this %C and the MAR from Eq. (4):

$$\begin{array}{}\text{(7)}& \text{CAR}\left(\mathrm{g}\mathrm{C}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\right)=\mathit{\%}{\mathrm{C}}_{\mathrm{core}}\cdot \text{MAR}\left(\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\right).\end{array}$$

Average CARs were calculated for the high marsh, the low marsh, and the entire region, using the eight cores with ²¹⁰Pb dating both for the whole cores as well as to the 30-year horizon (Table 2). Marsh-wide C accumulation rates (C accumulation) were calculated by multiplying CAR by the marsh area, for the high marsh, low marsh, and then for the full marsh (Table 2):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathrm{C}\text{accumulation}\left(\mathrm{g}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}\right)\\ & =\left[{\mathrm{CAR}}_{\mathrm{high}\mathrm{marsh}}\left(\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\right)\times {\mathrm{area}}_{\mathrm{high}\mathrm{marsh}}\left({\mathrm{m}}^{\mathrm{2}}\right)\right.\\ & \left.+{\mathrm{CAR}}_{\mathrm{low}\mathrm{marsh}}\left(\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\right)\times {\mathrm{area}}_{\mathrm{low}\mathrm{marsh}}\left({\mathrm{m}}^{\mathrm{2}}\right)\right].\end{array}\end{array}$$

Marsh-wide annual carbon accumulation was also estimated for the whole study area by multiplying the average CARs (high marsh, low marsh, and total) by the area of the seven marshes (Table 2).

2.8 Statistical analysis

We use a simple t test to test for significant differences in marsh characteristics (depth of peat base, carbon content, DBD, and SCDs), C stocks, and CARs between high and low marsh zones. We also compare base of peat C stocks, 20 cm C stocks, and 30-year C stocks, as well as for CARs estimated to 30 years and to the base of peat layer.

3 Results

3.1 Soil properties

With a few exceptions, marsh soils in Clayoquot Sound consisted of three layers separated by defined horizons of peat and sand/clay layers (Tables A1 and A2, Fig. A1 in the Appendix). In all cores, organic C concentrations were highest in the surface layers (10 %–45 %) and decreased to lowest values (∼0 %) in the deepest parts of the cores (Fig. 2).

Depths of refusal in the 34 cores ranged from 7 cm in a low marsh core of CBW to a maximum of 62 cm at CBE. The average (± standard error (SE)) corrected depth of the deepest peat layer overlying sand, gravel, or clay was 29 ± 3 cm for the high marsh sites and 17 ± 3 cm for the low marsh sites, with an overall average base of peat depth of 23 ± 2 cm (Fig. 3a). Differences in the deepest peat layer depths were statistically significant between the high and low marsh zones (p value <0.05). We note, however, that the 10 cores which did not penetrate into the sand, clay, or gravel layer were not used in these calculations.

While the average base of peat depth of high marsh cores is significantly higher than that of low marsh cores, no significant differences were found for average DBD, average %C, or average SCD (p>0.05) (Fig. 3b–d). Carbon contents averaged 21 ± 2 % in the high marsh and 18.5 ± 2 % in the low marsh cores and 20 ± 1 % overall (Fig. 3b) and did not vary significantly between the high and low marsh zones. DBDs and SCDs also showed similar average values and distributions in the high and low marsh cores (Fig. 3c and d). When corrected for compaction, soil carbon densities averaged 0.027 ± 0.002 g C cm⁻³ for all sites and ranged from 0.008 to 0.067 g C cm⁻³.

Figure 2 Percent carbon (vertical axis) versus depth (horizontal axis) for high marsh (black) and low marsh (grey) cores at all sites. Marsh name abbreviations are defined in Fig. 1. Depths have been corrected for core compaction.

Figure 3 Distributions of (a) deepest depth of peat overlying sand, gravel, or clay (cm); (b) average percent carbon; (c) average dry bulk density (corrected for compaction); and (d) average soil carbon density (corrected for compaction). Panels (b), (c), and (d) are estimated to the deepest peat depth for cores in the high marsh zones, low marsh zones, and all cores combined. Horizontal bar is the median; ends of the boxes represent the upper and lower quartiles; whiskers indicate minimum and maximum values, while dots indicate outliers. Differences between depths of deepest peat layer for high and low marshes are statistically significant (t=2.86, p value = 0.0086). Unpaired t tests show that the high and low marsh values are not statistically different for %C (t=1.03, p value = 0.31), DBD (t=0.25, p value = 0.80) , and SCD (t=0.91, p value = 0.37).

3.2 Carbon stocks and total carbon storage

The seven marshes ranged in size from 0.5 to 27 ha, with a total area of 47 ha (Table 1). The high marsh comprised 19 %–63 % of each individual marsh area and 58 % of the area of all seven marshes combined (Table 1).

The average C stock to the base peat layer for all cores was 67 ± 9 Mg C ha⁻¹ (mean ± SE). The base of peat C stock estimates for the high marshes (80 ± 14 Mg C ha⁻¹) were 1.5-fold higher than those found in low marsh cores (52 ± 8 Mg C ha⁻¹, unpaired t test, p<0.05) (Figs. 4 and 5, Table 1).

Similar results were obtained when C stocks were calculated only for the upper 20 cm (Fig. 4, Table 1). The average 20 cm C stock was 54 ± 5 Mg C ha⁻¹. High marshes averaged 64 ± 6 Mg C ha⁻¹, which was statistically higher than the low marsh average of 43 ± 7 Mg C ha⁻¹ (p<0.05). An average of 81 % of the total C stocks in a core accumulates in the top 20 cm of soil.

Summing C stocks to the 30-year horizon of the eight dated cores results in an average of 56 ± 14 Mg C ha⁻¹ (Table 1; Fig. 4). Using the 30-year horizon accentuates the statistically significant difference in C stocks accumulating in the high vs. low marsh zones (p<0.05), with high marsh C stocks (87 ± 16 Mg C ha⁻¹) up to ∼4 times higher than those found in low marshes (24 ± 2 Mg C ha⁻¹).

When combined with our estimates of marsh area, we calculate that the total C storage of these seven marshes was 4189 ± 303 Mg C, 69 % of which was stored in the high marsh (Table 2; Fig. 5c). For only the top 20 cm, this total is 2764 ± 149 Mg C, 63 % of which accumulated in the high marsh zones. Finally, over the last 30 years, 2852 ± 438 Mg C has accumulated in these marshes, 84 % of which accumulated in the high marsh zones (Table 1).

Figure 4 Carbon stock estimates for Clayoquot Sound for the high marsh zone (green), low marsh zone (blue), and combined (grey) for all cores. Estimates are made to the base of the peat layer, for the surface 20 cm, and for the most recent 30 years of accretion. Note that 30-year values are based on four cores each for high and low marshes (see Methods). Unpaired t tests show statistically significant differences between the high and low marsh C stocks calculated to the base of peat depth (t=2.91, p value = 0.0065), 20 cm C stocks (t=2.20, p value = 0.036) and 30 year C stocks (t=3.95, p value = 0.0074).

Figure 5 (a) High and low marsh carbon stocks (Mg C ha⁻¹) measured to the base of peat layer and for the last 30 years; (b) high and low marsh carbon accumulation rates (CARs, g C m⁻² yr⁻¹) measured to the base of peat layer and for the last 30 years; (c) total marsh carbon (Mg C) estimate to the base of peat layer for high and low marsh zones of each marsh; and (d) marsh-wide annual C accumulation (Mg C yr⁻¹) estimated to the base of peat layer and for the last 30 years. Paired differences between high and low marsh CARs are statistically significant for CARs calculated to the base of peat layer (t=4.283, p value = 0.0234) and for 30-year CARs (t=4.043, p value = 0.0272). (Statistical comparisons of C stocks are in Fig. 4).

Table 1 Carbon stocks.

Base of peat C stock was estimated by summing the mass of carbon from each 1 cm sample interval to the base of each core. Base of peat C stock values are averages of multiple cores, except where only one core was collected (e.g. high marsh KCS and low marsh GBK). The 30-year high and low marsh average C stocks are based only on CBE, CRF, GBK, and TMF. The 30-year marsh carbon (Mg C) is estimated by multiplying high marsh areas by average high marsh 30-year C stocks (2384 ± 438 Mg C) and low marsh area by average low marsh 30-year C stock (468 ± 19 Mg C), which is summed to get a total of 2852 ± 438 Mg C).

3.3 Age–depth relationships and mass accumulation rates over time

Activities of ²¹⁰Pb ranged from 2 to 372 Bq kg⁻¹, with maximum activities found at the surface. Activities of ²²⁶Ra ranged from 1.6–4.3 Bq kg⁻¹ and were subtracted from the ²¹⁰Pb_total to obtain the excess ²¹⁰Pb (²¹⁰Pb_xs). For most of the cores, the ²¹⁰Pb_sup level was not reached. Only GBK 1–4 and TMF 1–2 showed a complete ²¹⁰Pb_xs profile.

Activities of ²¹⁰Pb_xs showed an almost monotonic decline with depth suggesting that processes such as bioturbation, mixing, or disconformities have not disturbed the records and the logarithmic ²¹⁰Pb_xs profile plotted against the cumulative mass resulted in a linear profile (Fig. 6).

3.4 Mass accumulation rates and carbon accumulation rates

The “basal age”, defined as the age at the base of the deepest peat layer, ranged from 13 to 140 years (Table 2). At all sites, the basal ages are consistently older in low marsh cores (106 ± 17 years) when compared to high marsh cores (33 ± 9 years). As a result, the MARs were an order of magnitude higher in high marsh zones compared to low ones, with average MAR of 3978 ± 1231 and 399 ± 67 g m⁻² yr⁻¹ in the low marsh sediment, respectively, averaging 2188 ± 885 g m⁻² yr⁻¹ (range of 262–6977 g m⁻² yr⁻¹) for the whole marsh (Table 2).

The average base of peat CAR for the dated cores was 184 ± 50 g C m⁻² yr⁻¹ (range of 47 to 386 g C m⁻² yr⁻¹, Table 2). A similar value was obtained by using the 30-year average CARs, with an average of 178 ± 46 g C m⁻² yr⁻¹ (range of 59 to 394 g C m⁻² yr⁻¹) (Table 2; Fig. 5b). In the base of peat estimates, CARs in the four high marsh cores (average of 303 ± 45 g C m⁻² yr⁻¹) were significantly higher than the low marsh core average (63 ± 7 g C m⁻² yr⁻¹; p<0.05). This difference was also apparent when comparing the 30-year average CARs (280 ± 52 g C m⁻² yr⁻¹ and 76 ± 7 g C m⁻² yr⁻¹ for high and low marsh, respectively, p<0.05).

When we use the high and low marsh areas with the average of the CARs calculated to the base of the peat in each core, we find that the high marsh accumulates 83 ± 12 Mg C yr⁻¹ and the low marsh accumulates 13 ± 1.4 Mg C yr⁻¹ (Fig. 5c, Table 2). Combining these estimates produces an annual marsh-wide C accumulation of 96 ± 12 Mg C yr⁻¹, 86 % of which is accumulating annually in the high marsh zones. Similar results are obtained when using the 30-year CARs, with high marshes accumulating 85 % of marsh-wide annual C accumulation (77 ± 14 Mg C yr⁻¹) (Table 2).

Table 2 Core and marsh-wide carbon accumulation rates.

^a Marsh areas represent the sums for the high, low, and total marsh. ^b Marsh-wide carbon accumulation rates are estimated as the average (high, low, total) CAR multiplied by the marsh area (first column) for all seven marshes. ^c Associated DBDs, %C, and SCD values are in Table A1 in the Appendix. ^d Basal age reflects age determined in the deepest sediment layer containing peat overlying sand, clay, or gravel. ^e The total marsh-wide carbon accumulation values stated here are based on using the average total marsh CAR multiplied by the total marsh area. The values stated in the text are slightly higher because they are based on adding the low marsh and high marsh annual carbon accumulation values together. Please note that LSR denotes the linear sedimentation rate.

Figure 6 Natural log of excess ²¹⁰Pb activity versus cumulative mass from the surface. The equation in each graph shows the results from the regression analyses used for the calculation of the mean accumulation rate (g cm⁻² yr⁻¹) using the CF : CS model. Grey-shaded area represents 95 % confidence intervals.

4 Discussion

4.1 Carbon stocks

We calculate that the total C storage from salt marshes in Clayoquot Sound is 4189 ± 303 Mg C, with only 31 % accumulated in the low marsh (see Sect. 3.2). Similar, small contributions from the low marsh zone have been found in other estuarine systems (Connor et al., 2001; Adams et al., 2012). In nearby Boundary Bay, BC, Gailis et al. (2021) also found that the low marsh contributed a disproportionately small fraction (15 %) to the total marsh-wide C stocks relative to their area. As with Boundary Bay, we found that the small low marsh contribution was driven in part by the smaller proportion of low marsh area (42 % of total) in Clayoquot Sound. Gailis et al. (2021) also attributed lower C stocks in low marsh sediments to vegetation type and associated rooting depths (which are shallower in low marsh Salicornia spp.). Similar differences in vegetation are observed in the Clayoquot Sound marshes. Furthermore, when normalized to a constant depth horizon (20 cm), the C stocks in the Clayoquot Sound low marsh zones are 28 % lower than those in the high marsh zones (Fig. 4). However, the main drivers of the smaller C stocks are significantly shallower depths of the base of the peat layer in low marsh relative to the high marsh cores in Clayoquot Sound (Fig. 3a). High marsh peat layers are on average 71 % thicker than those found in the low marsh zones, suggesting that low marsh sediments are simply too shallow to accumulate as much C as the high marsh zones.

The C stocks in Clayoquot Sound marshes are less than one-third of the globally averaged estimate of 250 Mg C ha⁻¹ for salt marsh C stocks in the upper 50 cm (Chmura et al., 2003; Pendleton et al., 2012). This is true whether we consider the base of peat, 20 cm, or 30-year estimates (Table 1). One possible contributor to these discrepancies could be the use of different formulas relating %C to LOI. Chmura et al. (2003) utilized the formula of Craft et al. (1991) to convert LOI measurements to %C. In contrast, we used an empirical relationship based on measurements collected from Clayoquot Sound marsh samples. Comparison of these two regression equations suggest that estimates of %C are very similar for values of %C equal to or less than 30 %, but the %C estimates diverge for percentages above 30 % (Fig. A3), with calculated differences in %C exceeding 20 % at %LOIs above 80 %. However, we note that the interquartile (Q1 : Q3) range of our %C values fall between 4.4 % and 28.8 %, suggesting that most of our samples have %C values of less than 30 %, where the two equations produce similar results. To test the potential impact of the different equations, we compared C stocks (estimated to peat base) that were calculated using the two different %C–LOI relationships for the eight cores that were ²¹⁰Pb dated (Fig. A4). Using the Craft et al. (1991) regression inflates our C stock values by about 30 %, but this increase is not sufficiently large to account for the full difference between our C stocks and the global values in Chmura et al. (2003), which are 3 times greater than our estimated C stocks.

The C stocks in these seven Clayoquot Sound marshes are consistent with regional values found in Boundary Bay, BC (83 ± 30 and 39 ± 24 Mg C ha⁻¹ for high and low marsh, respectively, Gailis et al., 2021), in Stillaguamish, Washington (60 ± 4 Mg C ha⁻¹; Poppe and Rybczyk, 2021), and in Snohomish estuary, Washington (72 and 54 Mg C ha⁻¹ for the top 30 cm of two marshes, Crooks et al., 2014). However, these C stocks are an order of magnitude lower than estimates of total soil C mass in high and low marshes along the Pacific coast of the USA (543 ± 47 and 411 ± 70 Mg C ha⁻¹, respectively, Kauffman et al., 2020).

The lower overall C stocks in these Pacific marshes may in part be attributable to a lower overall C content: the average soil C density of the seven Clayoquot Sound marshes is 0.027 ± 0.002 g C cm⁻³, which is lower than regional estimates for the Oregon coast (0.034 ± 0.011 g C cm⁻³; Peck et al., 2020) and the global mean of 0.039 ± 0.003 g C cm⁻³ (Chmura et al., 2003). However, low overall C content is not likely to be the main cause of the lower C stocks overall because these SCDs are comparable to the means for the US (0.028 ± 0.013 g C cm⁻³; Holmquist et al., 2018) and Pacific Northwest marshes (0.023 ± 0.002 g C cm⁻³) in which C stocks are an order of magnitude higher than Clayoquot Sound (Kauffman et al., 2020).

A more likely cause is simply the shallow thickness of the Clayoquot Sound marshes, with the depth to the base of the deepest peat layer averaging 23 ± 2 cm (17 cm for low marsh and 29 for high marsh). Similarly shallow depths of refusal (penetrating into sand or clay layers) have been found in seagrass beds of Clayoquot Sound (Postlethwaite et al., 2018) as well as in Boundary Bay marsh (Gailis et al., 2021). In contrast, a sampling of seagrass and salt marsh systems along the US Pacific coasts suggest that these ecosystems have substantially thicker soil C depths, with >85 % having mean soil profile depths greater than 2 m (Kauffman et al., 2020). Because data represented in many syntheses come from regions with marshes that have depths of 1 m or greater (e.g. Macreadie et al., 2017), extrapolating the results to the total marsh area likely overestimates total C storage. Thus, we recommend that ongoing work to assess blue C contents globally begins to consider shallower depth horizons and/or age horizons as a means for creating a global comparison.

4.2 Carbon accumulation rates, comparisons of marsh zones and Pacific coast marshes

While marsh C stocks in Clayoquot Sound are substantially lower than the global and Pacific coast averages, the average CARs are actually higher than those found in other parts of the Pacific coast, especially when we separate measurements by marsh zone (Fig. 7). When all ²¹⁰Pb-dated sites are combined, we obtain an average CAR for the NE Pacific region of 112 ± 12 g C m⁻² yr⁻¹ (mean ± SE, n=18). While this value is comparable to the global average CAR of 91 ± 2.3 g C m⁻² yr⁻¹ (Kennedy et al., 2013), both are substantially lower than the regionally averaged CAR of 184 ± 50 g C m⁻² yr⁻¹ for Clayoquot Sound (Fig. 7c). When we consider only high marsh cores, the average CAR for Clayoquot Sound further increases to 303 ± 45 g C m⁻² yr⁻¹, substantially larger than the NE Pacific high marsh average of 145 ± 26 g C m⁻² yr⁻¹ (Fig. 7a). In fact, the NE Pacific high marsh average is driven substantially by the inclusion of Clayoquot Sound marshes and drops to 96 ± 12 g C m⁻² yr⁻¹ when Clayoquot Sound marshes are excluded. This comparison suggests that while the overall C storage in Clayoquot Sound marshes is low, the current rates of accumulation are much higher than the regional average, particularly in the high marsh zones. Thus, assuming these rates of C accumulation continue, high marsh zones of Clayoquot Sound have the potential to play a more important role as C sinks in the future.

One effort has been made to assemble a global compilation and synthesis of CARs within salt marsh ecosystems (Ouyang and Lee, 2014). Here we compare our results with this compilation with the caveats that (a) Ouyang and Lee (2014) relied heavily on ¹³⁷Cs dating or marker horizon methods for their estimates and (b) there appear to be some instances of double counting of sites (e.g. averages of three sites were included as a fourth site) and minor issues with attribution of some geographic locations. That said, we note that our new average estimate for NE Pacific CARs of 112 ± 12 g C m⁻² yr⁻¹ is substantially lower than estimates compiled by Ouyang and Lee (2014), whether considering the global average (245 g C m⁻² yr⁻¹), the NE Pacific average (167 g C m⁻² yr⁻¹; n=6), the NW Atlantic region (172 g C m⁻² yr⁻¹; n=64; 35.0–47.4^∘ N), or the subset of cores from Atlantic Canada (188 g C m⁻² yr⁻¹; n=40; 43.6–47.4^∘ N). For the NE Pacific region, this compilation predominantly included sites where ¹³⁷Cs dating was used to estimate CARs (Ouyang and Lee, 2014). The remaining NE Pacific CAR studies (n=2) used marker horizon methods which, when included, increase the average CAR to 174 ± 128 g C m⁻² yr⁻¹. Marker horizon and ¹³⁷Cs dating methods are shown to produce CARs higher than those calculated from ²¹⁰Pb (Chmura et al., 2003; Callaway et al., 2012). This is due to gradual decomposition in the sediment column, averaging CAR over a shorter time period, and to a shallower depth, oversampling high C soils vis-à-vis lower C soils at depth. As most of the CAR data points in the NW Atlantic region in Ouyang and Lee's 2014 summary were also calculated using ¹³⁷Cs dating or marker horizon methods (n=62), these also represent roughly equivalent values to Clayoquot Sound if the observed bias is considered.

Differences between Clayoquot Sound and regionally averaged CARs are not apparent for the low marsh zone, for which the Clayoquot Sound average of 63 ± 7 g C m⁻² yr⁻¹ is statistically similar to the regional estimate of 74 ± 6 g C m⁻² yr⁻¹ (Fig. 7b). Within the context of Clayoquot Sound, the low marsh zone exhibits substantially lower C stocks and accumulation rates than the high marsh zones. Basal ages (age of the deepest peat layer overlying sand) are substantially older in low marsh sediments when compared to the high marsh, resulting in significantly lower sedimentation rates and CARs.

Our observation of significantly lower CARs in the low marsh compared to the high marsh could result from several factors. As mentioned previously, the vegetation in the low marsh zone is dominated by Salicornia spp., whose lower productivity and shallow root systems could be producing and trapping less C (Kelleway et al., 2017; Gailis et al., 2021). A second potential contributor to slow accumulation rates of C could be the relatively low sediment discharge from terrestrial sources such as the Kennedy River (Nuwer and Kiel, 2005). Previous work in Clayoquot Sound attributed low carbon accumulation rates in eelgrass, in part, to relatively clear waters (Secchi disc measurements) and $\mathrm{C}/\mathrm{N}$ ratios that reflect marine rather than terrestrial sediment contributions (Postlethwaite et al., 2018).

A third important feature of the Clayoquot Sound marshes is that the region is experiencing vertical land motion that counteracts the effects of sea level rise (Mazzotti et al., 2007, 2008; Montillet et al., 2018). Previous work suggests that, in many cases, marshes are responding to rising sea levels by migrating inland (Kirwan et al., 2016). In contrast, in our study region the nearby town of Tofino experiences uplift of 2.86 mm yr⁻¹, such that relative sea level is dropping by 1.08–1.15 mm yr⁻¹ (Montillet et al., 2018). In fact, sea level change along the Pacific coast of Canada where these marshes are developing averages −0.76 ± 1.32 mm yr⁻¹ (James et al., 2014), suggesting that, while some areas along the Pacific Canadian coast are exposed to sea level rise, average conditions for regional salt marshes are likely comparable to what is being experienced in Clayoquot Sound. Recent work along the Oregon coast suggests that relative sea level rise is a dominant control on vertical accretion rates in tidal marsh systems, and areas experiencing negative sea level rise (or land uplift) have the lowest rates of C accumulation (Peck et al., 2020). The negative sea level rise and relatively low C accumulation rates, particularly in the low marsh settings of Clayoquot Sound and Pacific Rim National Park, are consistent with this result.

A final discriminating characteristic between the high and low marsh cores in Clayoquot Sound are the ²¹⁰Pb inventories (Fig. A5), which suggest that the high and low marsh zones are subjected to different erosional processes. Given that these cores all come from the same area one would assume that all the cores should have similar inputs of ²¹⁰Pb (Olid et al., 2010; Arias-Ortiz et al., 2018). However, a comparison of the ²¹⁰Pb inventories shows that the high marshes accumulate much more ²¹⁰Pb than the low marshes. This difference, together with the lower concentrations of the ²¹⁰Pb in the upper layers of the low marsh cores, could result from focusing/soil formation causing greater ²¹⁰Pb accumulation in the high marshes and/or reduced accumulation in the low marshes due to erosion and winnowing processes. One possible explanation here is that the extensive, shallow flats in the low marsh regions of Clayoquot Sound are subjected to more intensive tidal erosion and therefore removal of soil organic material, which results in lower accumulation rates relative to the high marsh zones. Further studies are needed to estimate the significance of these drivers within the region.

Figure 7 Comparison of Clayoquot Sound CARs with other salt marshes with ²¹⁰Pb dating, for whole-core (black circles) and 30-year (red triangle) estimates. Individual Clayoquot Sound marshes are indicated with open symbols followed by a regional mean (± SE). Additional site data taken from (1) Boundary Bay, BC, Canada (Gailis et al., 2021); (2) Stillaguamish, WA, USA (Poppe and Rybczyk, 2021); (3) Snohomish Estuary, WA, USA (Crooks et al., 2014); (4) Animal Slough (AS) in southern Puget Sound, WA, USA (note that red triangle is 50-year average; Drexler et al., 2019); (5) coastal OR, USA (Peck et al., 2020); and (6) San Francisco Bay, CA, USA (Callaway et al., 2012). Global and other regional averages are taken from (7) Kennedy et al. (2013) and (8) Ouyang and Lee (2014); only a fraction of data (<10 %) in (8) were dated using ²¹⁰Pb. Horizontal line represents Kennedy et al. (2013) global average.

4.3 Comparison of whole-core and 30-year CARs

Estimating CARs to a consistent age horizon of 30 years allowed us to compare C accumulation between marshes and marsh zones, but it also reveals interesting information about differences in CAR depending on the approach used to estimate these values. Previous comparisons in nearby Boundary Bay marsh showed that CARs were 20 % higher during the last 30 years when compared to the whole-core estimates (Gailis et al., 2021). Because %C decreases with depth in all 34 cores from Clayoquot Sound, we might anticipate that CARs for the most recent years will always result in a higher CAR than a measurement of C stock to the depth of minimum ²¹⁰Pb detectable activity. This can be partially attributed to remineralization of soil C over time, which would decrease overall %C in the deeper portions of the sediment cores (Callaway et al., 2012; Johannessen and MacDonald 2016). In Clayoquot Sound, however, the CARs estimated for the last 30 years were not substantially different from our estimates made to the base of the peat layer (Table 2), particularly in the high marsh zone where the average age for the base of the peat layer was 33 ± 9 years. Slightly larger differences were seen in the slowly accumulating low marsh zones, where the 30-year CAR was about 15 % higher. Overall, the young ages of these marshes result in small differences between the two estimates.

The examination of 30-year CARs as well as CARs in relatively young marshes could have interesting implications for future blue C assessments. Many marsh restoration projects in which C is being measured are currently on timescales of 10 to 30 years past the initiation of restoration (e.g. Abbott et al., 2019; Drexler et al., 2019). Vegetational succession and C storage in a restored tidal marsh can take over 10 years (Havens et al., 2002; Craft et al., 2003), and even a 30-year monitoring period would capture a period of gradually improving C sequestration potential in restored tidal marsh soils (e.g. Howe et al., 2009). However, the inflation of CARs in younger sediments in the Clayoquot Sound and Boundary Bay sediments could provide a caution that CARs measured from relatively young marker horizons are at risk of overestimating total C in sediments, especially if these CARs are used to project C sequestration over longer time periods. Several recent studies have shown substantially higher CARs in restored vs. natural marshes after approximately 10 years (e.g. Poppe and Rybczyk, 2021; Drexler et al., 2019). The causes of these inflated CARs can be very site-specific, depending on the types of colonizing vegetation (Kelleway et al., 2017) and the relative contributions of allochthonous mineral materials (Drexler et al., 2019). In the case of Clayoquot Sound, it is possible that these high CARs will be sustained, making this marsh extremely important as a regional C sink. However, we also cannot rule out the possibility that longer-term processes, such as remineralization and sediment redistribution, could reduce the actual amount of C stored as the marshes age, such that the accumulation rates become more comparable to those found on other parts of the Pacific coast (Fig. 7).

4.4 Comparison with carbon storage in boreal forest ecosystems

The relevance of blue C storage potential to climate change mitigation depends on the scale over which it is considered. In the literature, the carbon accumulation potential of blue carbon is often compared with terrestrial ecosystems, and it is often stated that, per unit area, blue C ecosystems accumulate carbon at higher estimated rates than terrestrial forest soils (McLeod et al., 2011).

A first-order comparison for Canadian ecosystems shows the same pattern of higher carbon uptake rates per unit area in tidal wetlands when compared with measurements of net ecosystem productivity of terrestrial forests. Forested areas in Canada have been estimated to take up carbon at rates ranging from 35 g C m⁻² yr⁻¹ (Canada-wide estimate, Stinson et al., 2011) to 63 g C m⁻² yr⁻¹ (British Columbia, Peng et al., 2014). The 30-year CAR for salt marshes in Clayoquot Sound is approximately 2–7 times higher, averaging 178 ± 46 g C m⁻² yr⁻¹. Note that we considered forests' net ecosystem productivity, which is higher than C uptake rates of 4.6 ± 2.1 g C m⁻² yr⁻¹ for boreal forest soils cited by McLeod et al. (2011).

Consideration of the areal extent of each of these ecosystems is extremely important for contextualizing blue C uptake among Canadian ecosystems. Estimates of total salt marsh area in Canada range from 44 000 ha (Bridgham et al., 2006) to 111 274 ha (Mcowen et al., 2017). In contrast, the boreal forest ecosystem is one of Canada's largest terrestrial biomes and encompasses approximately 270×10⁶ ha (Kurz et al., 2013). Our Clayoquot Sound data represent only a small area of a single region of the west coast, but if we assume our CAR estimate of 184 ± 50 g C m⁻² yr⁻¹ from Clayoquot Sound approximates the average for all tidal salt marshes in Canada, Canada's marshes would accumulate somewhere between 81 000 and 205 000 Mg C yr⁻¹. This annual carbon accumulation from Canadian salt marshes is between 0.4 % and 2.9 % of the carbon estimated to accumulate annually by Canadian boreal forests (6 750 000–18 090 000 Mg C yr⁻¹).

Our back-of-the-envelope comparison is not intended to provide a formal estimate of carbon sequestration in Canadian salt marshes but instead highlights several important assumptions and uncertainties that are important to large-scale assessments of blue carbon sequestration. The first is that regional-to-global estimates of blue C accumulation rates depend on regionally specific measurements of salt marsh C stocks and accumulation rates, and extrapolation of blue C accumulation rates based on globally derived averages (or averages from regions other than the one in question) can lead to large uncertainties in blue C accumulation rates. Our calculation above assumes that the CAR estimates from Clayoquot Sound are representative of the pan-Canadian average for CARs in salt marshes. Using the published value for Atlantic Canada (188 g C m⁻² yr⁻¹, Ouyang and Lee, 2014) would not substantially change our calculation. However, if the average pan-Canadian CAR is closer to the global estimate of 91 ± 2.3 g C m⁻² yr⁻¹ (Kennedy et al., 2013), our estimate of annual carbon accumulation would halve. Thus, a pan-Canadian estimate of CAR (and its associated uncertainty) would be useful for improving the estimated carbon sequestration rates in Canadian salt marshes.

Second, our comparison with C sequestration in Canadian forests is based on net ecosystem productivity estimates from 2011 and 2013. With climate change, the increased prevalence of natural disturbances such as forest fires and insect outbreaks is projected to result in episodic losses of carbon from the above-ground and soil carbon pools (e.g. Stinson et al., 2011; Walker et al., 2019). Between 1990–2008, 24 ± 19 Tg C yr⁻¹ was removed from the forest carbon pool by fire (Kurz et al., 2013), with British Columbia forests becoming a sustained source of CO₂ since 2003 (e.g. Giles-Hansen and Wei, 2022). The degree to which climate change will impact forest carbon sequestration in future will depend on several factors, including increases in disturbance frequency and intensity (Walker et al., 2019) and subsequent recovery and shifts in forest species dominance (e.g. Mack et al., 2021), as well as management practices (Price et al., 2013). At the same time, as sea level rise is currently a minimal threat to the Canadian Pacific salt marshes, they are likely to continue to function as efficient C sinks despite global warming.

Finally, in addition to the regionally specific nature of processes controlling blue C accumulation, we note that estimation of the areal extent of these marshes varies by an order of magnitude within Canada alone, and determinations of the relative proportions and high-to-low marsh settings have not been quantified. We note, for example, that while large-scale databases suggest that the areal extent of salt marsh in Boundary Bay, BC, is 1207 ha (CEC, 2016a), our field research in Boundary Bay has determined that the actual marsh area is only one-quarter the size at 275 ha (Gailis et al., 2021). This suggests that estimation of blue C potential in policy-specific contexts requires a more detailed determination of the areal extent of these marshes, especially when considering provincial-to-national scales. We further suggest that the small area over which blue C accumulates likely influences the scale over which it is best considered for natural land carbon sequestration and accounting. Blue C accumulation occurs in only 0.016 %–0.1 % of the land area of boreal forests nationwide. On a global, national, and even provincial scale, total blue C sequestration in salt marshes is relatively small compared to the large areal expanses of forest. However, the very high, per-unit area sequestration of carbon in salt marsh ecosystems can increase its importance at local-to-regional scales and should be considered an important co-benefit to restoration activities and conservation for other purposes.

4.5 Future work

While this study has provided a useful first step in quantifying blue C stocks, storage, and accumulation rates, we note the following considerations in estimating the potential of blue C within this region. First, our study has estimated blue C in 47 ha of high and low marsh areas of Clayoquot Sound. We expect that these mesotidal estuarine marshes, often constrained in size by surrounding topography, are typical of the marshes found on the Pacific coast of British Columbia, the total estimated area of which is 60 km², or 6000 ha (Ryder et al., 2007). However, given variability in riverine discharge and tectonic settings along the coast, we anticipate that additional measurements (additional marshes, with differing riverine and sedimentary inputs, and variable uplift setting) will be necessary to understand the variability in carbon accumulation along the Pacific coast of Canada. Second, our study has focused on estimation of carbon stocks and accumulation rates within these marshes, and full estimation of blue C potential will ultimately require estimation of greenhouse gas fluxes of carbon dioxide, methane, and nitrous oxide within these marshes (e.g. Chmura et al., 2011). The brackish nature of some of these marshes could mean that some emit substantial contributions of methane, which may counter their effectiveness as C sinks (e.g. Poffenbarger et al., 2011; Abdul-Aziz et al., 2018; Huang et al., 2019; Huertas et al., 2019; Li et al., 2021; Wei et al., 2020). Finally, this study has revealed the importance of classifying high vs. low marsh zones to best quantify carbon stocks and accumulation rates, as the different vegetation growth and sedimentary processes operating in these zones can critically influence CARs and therefore estimates of C sequestration. Future work thus needs to involve better quantification of the areal extent and classification of salt marshes, within Clayoquot Sound, along the Pacific coast, and in Canada as a whole.

5 Conclusions

Our work helps address the data gap for North American blue C by producing CAR estimates using ²¹⁰Pb-derived mass accumulation rates of salt marsh sediments, showing CARs comparable to other salt marsh locations on the northeast Pacific coast. Importantly, this work demonstrates that the marshes of Clayoquot Sound along the Pacific coast of Canada exhibit substantially smaller C stocks when compared with other regional salt marshes to the south, along the Pacific coast of the United States. We attribute these lower stocks to substantially shallower (<50 cm) depths of marsh initiation in all high and low marsh settings. At the same time, while carbon stocks are lower than regional and global averages, the rates of carbon accumulation in these marshes are comparable and even higher than regional estimates. We show that both 30-year CARs and CARs estimated to the peat base are substantially greater than regional averages generated using only ²¹⁰Pb-derived mass accumulation rates. However, significant variability exists between marshes, and different radioisotope dating methods in CAR measurement create uncertainty. Our results support existing findings that blue C soils accumulate carbon at high rates, limited only by the spatial extent of such ecosystems. These results provide a natural analogue for monitoring carbon accretion over timescales relevant for restoration activities.

Appendix A: Supplemental tables and figures

Table A1 Core locations, soil carbon characteristics, and carbon stocks.

^a Avg DBD values were not corrected for compaction (see Table A3 for corrected DBDs), but SCD values have been corrected for compaction. ^b Cores that do not extend to 20 cm and do not penetrate a sand/gravel/clay layer at the base were excluded from the 20 cm calculation. ^c GBK cores were extracted using a steel sledge corer that penetrated 15–30 cm into the sand and clay layers. The total core lengths were 46 (88), 60 (100), 59 (66), and 24 (37) cm for GBK 1–1, 1–2, 1–3, and 1–4, respectively. Depths provided here are just below the transition from peat to sand, 2 cm into the minimal %C layer.

Table A2 Downcore distribution of ²¹⁰Pb and ²²⁶Ra in eight cores used for ²¹⁰Pb dating.

^a Measurements from Core Scientific International. ^b Measurements from Flett Research Ltd. ^c Measurements from MyCore Scientific Laboratory. ^d “Soil” represents a black, fluffy surface organic layer overlying peat in several cores. Brown peat layer below contained extensive, decaying plant roots. The transition between these two layers was visually distinct and also reflected in changes in dry bulk density. ^e Depth of 30-year horizon. ^f Depth basal age (age determined in deepest sediment layer containing peat overlying sand, clay, or gravel). Note that the 30-year horizon can occur below the depth of the basal age. Bold values represent the depth levels of the 30-year and basal-peat depth horizons, used for calculating C stocks.

Figure A1 Core stratigraphies. Dated high marsh cores (asterisks) were CBE 1–1, CRF 2–1, GBK 1–2, and TMF 2–1; dated low marsh cores were CBE 1–5, CRF 3–2, GBK 1–4, and TMF 1–2. Corrected depth (cm) reflects correction for core compaction. Surface organic layer (“soil” in table ) overlying peat in several cores was black, fluffy, highly organic material overlying the brown peat layer below, which contained extensive plant roots that had not entirely decayed yet. The transition between these two layers was visually very distinct and also reflected by changes in dry bulk density.

Figure A2 %LOI-to-%C relationship. Relationship between measured LOI values and calculated %C, using elemental analyzer (EA) data on set of 93 subsamples. Measurements from core CBE 1–5 were not used for calculating this relationship due to suspected measurement error.

Figure A3 Comparison of %C calculation using the empirical regression determined for our study area (blue) and the model of Craft et al. (1991) (orange).

Figure A4 Comparison of C stocks (Mg C ha⁻¹) estimated using the Craft et al. (1991) equation versus C stocks estimated using this study's empirical relationship for southern BC, for eight cores from our study region. Comparison suggests that using the Craft et al. (1991) relationship would produce C stocks that are ∼32 % greater than our estimates. This difference, while substantial, would not account for global C stocks that are 3 times larger than those found in Clayoquot Sound salt marshes.

Figure A5 ²¹⁰Pb inventory (Bq m⁻²) in high and low marsh cores.

Appendix B: Ground-truthing marsh zone designations: detrended correspondence analysis

We used colour orthophotos to visually determine the areas of marsh zones using ArcMap 10.5.1 area measurement tools. A vegetation survey also determined the marsh zone of each core sample site, and vegetation data were examined using a detrended correspondence analysis to verify that vegetation assemblages reflected distinct low and high marsh zones identified from the orthophotos (Hill and Gauch, 1980).

Our comparison showed that differentiating high and low marsh using orthophotos matched with field vegetation data for 32 of 34 (94 %) cores. Both CRF 1–2 and CRF 2–2 were classified as low marsh by vegetation survey but fell within the high marsh using the visual orthophotography method. These cores lie 16 m (CRF 1–2) and 12 m (CRF 2–2) away from the low marsh boundary as measured using orthophotos, which is less than their distances from the nearest high marsh cores (17 and 23 m, respectively). All other cores fell within the correct marsh zone.

We also conducted a detrended correspondence analysis (Hill and Gauch, 1980) of vegetation data to determine that vegetation encountered in the vicinity of the core samples reflected distinct marsh zones using Canoco v4.5 software. We classified marsh zones by presence/absence and percent cover of low marsh or high marsh vegetation (e.g. Deur, 2000). The analysis showed an accurate fit of low marsh cores with low marsh vegetation and high marsh cores with high marsh vegetation, plus the addition of a somewhat indistinct, third cluster of vegetation possibly representing the upland vegetation. This square-root-transformed model accounted for 33.2 % of all variance in the vegetation dataset (sum of eigenvalues = 3.29). Cores with low marsh vegetation clustered together, while high marsh cores clustered separately. An additional, slightly distinct third cluster of upland vegetation indicates that some high marsh cores may have been extracted from near the boundary with a freshwater-dominated upland vegetation or salt-tolerant meadow. The distinction between a salt marsh and a bordering freshwater area has complicated efforts to classify marshes by salinity (Duarte et al., 2013), but this result shows that clustering of vegetation type corresponds reasonably well with each site's designation as high or low marsh.

Figure B1 Detrended correspondence analysis results for marsh vegetation data. Low marsh cores (top, purple squares) corresponded reasonably well with vegetation identified as low marsh and high marsh cores corresponded with a distinct cluster of high marsh vegetation. The far bottom right may indicate a population of less salt-tolerant, upland vegetation, but it is indistinct from the high marsh.

Figure B2 Reference used for determining marsh stratum based on vegetation, adapted from Deur (2000). Dashed vertical lines from left to right represent mean time (1 m), mean lower high water (1.5 m), mean higher high water (1.85 m), and maximum water level (2.85 m).