Spatial and temporal variability in blue carbon accumulation in the largest salt marsh in British Columbia, Canada

Preserving blue carbon ecosystems, such as salt marshes, for climate change mitigation requires quantifying their carbon (C) dynamics. Boundary Bay (BB) marsh is a 222-ha salt marsh in southern British Columbia, Canada, where construction began in 2023 to enhance the marsh as a natural defense against coastal flooding. This study provides a baseline understanding of C storage and sequestration in Boundary Bay marsh prior to foreshore enhancement. We collected 18 sediment cores and vegetation surveys across the middle (BBM), eastern (BBE), and Mud Bay (MB) areas of the marsh, along with 128 depth profiles (i.e., field measurements of marsh thickness to refusal) from BBM. We combined C measurements with ²¹⁰Pb chronologies, in addition to existing data from western Boundary Bay (BBW), to estimate C stocks (g C m^-2) and accumulation rates (g C m^-2 yr^-1) for the entire marsh. Total C stocks averaged 71 ± 37 Mg C ha^-1 for high marsh and 41 ± 36 Mg C ha^-1 for low marsh, with higher values in western Boundary Bay (BBW, BBM) compared to the east (BBE, MB). Total C storage (Mg C) at Boundary Bay was 17,360 ± 4,960 Mg C, with the western marsh (BBW, BBM) comprising 84% of the total. The C accumulation rates (CAR) for Boundary Bay marsh averaged 80 ± 45 g C m^-2 yr^-1, comparable to regional averages on the Pacific coast of North America. However, large spatial variability exists, with significantly lower average CARs in the east (35 ± 11 g C m^-2 yr^-1). Historical aerial photographs indicate that the eastern marsh area (BBE, MB) decreased by ~35% while BBW expanded by ~20% since 1930. These contrasting trends suggest dynamism in marsh development, likely driven by environmental factors and human influence. This work highlights the high spatial and temporal dynamics of blue C ecosystems, especially in urban settings, and how decadal changes induced by human activities could influence their short-term (years to decades) C storage capacity, with potential consequences for long-term (centuries to millennial) C sequestration.

1 Introduction

“Blue carbon” refers to the organic carbon (C) captured and stored by vegetated coastal ecosystems (such as tidal marshes, mangroves, and seagrasses) and in the ocean. Among these systems, tidal salt marshes are particularly important for their capacity to sequester and store C over long timescales (Chmura et al., 2003; Mcleod et al., 2011; Duarte et al., 2013; Howard et al., 2014; Macreadie et al., 2021), despite covering only ~0.04% of global land area between 60°S and 60°N (Worthington et al., 2024). Due to their capacity to mitigate climate change by sequestering C, restoration and conservation of tidal salt marshes have received increasing attention as a natural method to compensate for anthropogenic global C emissions (Janousek et al., 2021; Irving et al., 2011). Salt marshes provide multiple important ecosystem services (e.g. filtration, storm surge protection, biodiversity, C sequestration), which make them potentially important contributors to natural climate solutions (NCS). NCSs include a suite of improved management and restoration actions that increase C storage or reduce greenhouse gas emissions in landscapes and wetlands across the globe, providing additional climate mitigation beyond business as usual (Drever et al., 2021; ECCC, 2021). However, the inclusion of salt marshes – coastal wetlands rich in C - in NCSs is limited due to the lack of information on spatial and temporal estimates of salt marsh coverage and C storage. Boundary Bay salt marsh represents the largest salt marsh in British Columbia, and eastern Boundary Bay currently lacks specific knowledge about marsh area, C stocks, C accumulation rates (CARs), and historic processes of marsh expansion or erosion, all of which could influence the marsh’s capacity to serve as a long-term C sink.

Boundary Bay salt marsh, which spans more than 200 ha, is an excellent example of a coastal tidal marsh ecosystem that is being considered for application of an NCS. The bay provides critical habitat for many aquatic organisms and is considered a potentially important blue carbon ecosystem (Gailis et al., 2021; Dashtgard, 2011). Environment and Climate Change Canada has committed to help fund a “living dike” pilot project in Delta, Surrey, and on the traditional territories of the Semiahmoo First Nation (City of Surrey, 2021; Readshaw et al., 2018). A living dike is a type of flood protection structure that uses a naturally existing salt marsh to enhance coastal protection against storm surges and flooding.

The C sequestration potential of Boundary Bay salt marsh has not been considered through the implementation of the living dike project. While previous work has quantified the areal extent, C stocks, and CARs in the far western portion of Boundary Bay (Gailis et al., 2021), data on C dynamics are still lacking for the eastern portion of the marsh. Filling in this information gap is an important next step towards a more holistic understanding of temporal variation in the areal extent of the Boundary Bay marsh ecosystem and the different processes controlling C storage in the high and low marsh zones.

This study aims to map the total areal extent, volume, C stocks, and CARs in Boundary Bay marsh. Specifically, we address the following research questions: i) What is the current areal extent and volume of Boundary Bay marsh (and how does this area compare with historical estimates of marsh extent)? ii) What are the total C stocks and how do they vary across different subregions of the marsh? iii) What are the rates of carbon accumulation (CARs), and how do these compare with previous studies? By framing the study in this way, we aim to provide a more comprehensive quantification of C dynamics in Boundary Bay and to contribute to broader efforts of evaluating blue carbon ecosystems.

2 Materials and methods

2.1 Study area

Boundary Bay is a shallow, sheltered bay situated on the Canada-United States border in southwestern British Columbia (BC), within the traditional territory of the Semiahmoo Nation and adjacent to the municipalities of Delta to its west and Surrey to its east. Prior to the late 1950s, the shoreline of Boundary Bay was one of the most productive shellfish harvesting locations in the Pacific coast, mainly used for oyster farming. In early 1960’s, operation licences were suspended due to sewage pollution from the Serpentine and Nicomekl rivers (Readshaw et al., 2018). Boundary Bay has been a designated wildlife management area since 1995 (Government of British Columbia, 2020).

To its north, the entire length of the urban marsh is bound by a 26 km man-made dike built in 1895 and renovated in 1948 following the Fraser Valley Flood (Shepperd, 1981). The dike is adjacent to agricultural and residential land, greenhouses, a golf course, a domestic airport, and a highway. The Boundary Bay salt marsh spans 19 km from western Boundary Bay to Mud Bay and extends approximately 0.5 km at its widest points, measured perpendicular to the dike in its western and eastern portions.

Boundary Bay receives freshwater and sediment inflows from the Serpentine and Nicomekl Rivers, which lie to the northeast of Mud Bay. These small rivers originate ~35 kilometers inland and discharge a minor amount of sediment into the bay (Swinbanks and Murray, 1981). To the north, the Fraser River discharges large amounts of freshwater and sediment into the Strait of Georgia. However, freshwater and sediment inputs from the Fraser River into Boundary Bay are limited by the barrier of Point Roberts peninsula on the west side of the bay. The Point Roberts Peninsula also protects Boundary Bay from westerly and southwesterly winds. However, erosion from the Pleistocene cliffs at Point Roberts is a present-day source of sediments to the Bay, which include gravel, sand, silty sand, mud, peat, and driftwood/shell accumulations (Shepperd, 1981; Swinbanks and Murray, 1981). The circulation in the bay is affected by prevailing, large-wave-generating winds from the east and south-east, while tidal currents flow in from the south-west (Dashtgard, 2011). The bedrock beneath the bay is known to be Pliocene/Miocene formations (Shepperd, 1981).

The rates of sedimentation at Boundary Bay are low (0.42 mm yr^-1) with relatively clear, saline (24–29 ppt) waters (Swinbanks and Murray, 1981). The western portion of the marsh is expanding (Kellerhalls and Murray, 1969; Gailis et al., 2021) while most of the mid and eastern portions are thought to be receding due to the water current and wind dynamics (Swinbanks and Murray, 1981). The tides in the Bay are a mixed semidiurnal type, with two high and two low tides every day (Shepperd, 1981). The mean tidal range is 2.7 m, with a maximum spring tide of 4.1 m and minimum neap tide of 1.5 m (Swinbanks and Murray, 1981).

For the purposes of this study, the Bay was divided into four study areas: western, middle, and eastern Boundary Bay (BBW, BBM, BBE, respectively) as well as Mud Bay (MB), all located within the cities of Surrey and Delta, BC (Figure 1). First, BBW is in the northwestern tip of Boundary Bay, to the east of 64^th St, Delta, BC, that was previously sampled in 2014-2018 (Gailis et al., 2021). The salt marsh at this study site is 6 km long at its leading edge and 0.5 km wide at its widest point measured perpendicular to the dike. Second, BBM is located south of 72^nd St in Delta, BC, just to the east of BBW. Here, the salt marsh is 6 km long at its leading edge and 0.4 km wide at its widest point measured perpendicular to the dike. To the east, there is no marsh present for 3.2 km in the center portion of the bay due to water currents and wind patterns (Dashtgard, 2011). Third, BBE represents the second largest portion of the study area and is located on the eastern side of the bay in Delta, BC, to the south of the Highway 99 loop that was built in 1942. The salt marsh at this study site is 5 km long at its leading edge and 0.5 km wide at its widest point measured perpendicular to the dike. Finally, MB is in the far, northeast corner of Boundary Bay, in Mud Bay Park in Surrey, BC. The salt marsh at this study area is 2.0 km long at its leading edge and 0.1 km wide at its widest point measured perpendicular to the dike. MB is close to the freshwater influx from the Serpentine and Nicomekl Rivers. The marsh in this area is underdeveloped due to the low sediment input from the rivers and is thought to be receding due to eroding water currents (Dashtgard, 2011). A high marsh area was not present due to the small size of the marsh.

Figure 1

Figure 1. Boundary Bay marsh in the municipalities of Delta and Surrey, British Columbia, Canada. Map inserts show locations of four study areas of western (BBW), middle (BBM), eastern (BBE) and Mud Bay (MB) portions of Boundary Bay. Dated cores are indicated by solid red (high marsh) and solid blue (low marsh) dots. White dots indicate all other cores. Year-round pore water salinity measurements were collected at sites marked by white asterisks. Base map source: ArcGIS 2021.

2.2 Marsh area and volume

2.2.1 Marsh area

We estimated the area of high and low marsh zones using Google Earth Pro 2021 tools with a Google satellite base map compiled from images taken on June 12, 2019 (Google Earth Pro, 2019). High and low marsh zones were delineated by eye using differences in vegetation color on Google Satellite base map imagery (Figure 1). Low-marsh plant assemblages were a lighter shade of green compared to the darker green of the high marsh plant assemblages (Chastain et al., 2022; Gailis et al., 2021). Gailis et al. (2021) determined that using vegetation surveys aligned perfectly with the visual delineations made using satellite color. This delineation was further verified by vegetation surveys at the 18 sediment coring sites during the 2020–2021 sampling period. Therefore, all areas not containing sediment cores were delineated only by eye using 2019 Google Satellite base map vegetation coloration (Google Earth Pro, 2019).

High and low marsh areas for BBW, BBM, BBE and MB were calculated separately and then summed up to determine the total marsh area for Boundary Bay. The 140-ha study area previously identified by Gailis et al. (2021) as the “western marsh” was divided into BBW (80 ha) and BBM (84 ha) sites for this study. Additionally, a 24-ha low marsh area adjacent to the eastern edge of the marsh studied by Gailis et al. (2021) was added to our BBM study area to more accurately represent the total marsh extent.

2.2.2 Marsh volume

To estimate the total volume of the marsh, marsh depths were measured with their associated GPS coordinates at 128 sampling sites throughout BBE (the so-called “depth profiles”, Supplementary Table S3). Twenty transects separated every 50 m were selected; the seven sampling sites along each transect were approximately 20 m apart. The depth measurements were collected by inserting a 6-ft plasticized metal garden stick down the marsh until it hit the depth of refusal. Depth of refusal is a reasonable estimate of organic soil thickness because it assumes that organic soil is easier to penetrate than underlying sands and/or bedrock (Howard et al., 2014). For the purposes of this study, we assume that the underlying sand layer was formed prior to marsh initiation (Howard et al., 2014; Gailis et al., 2021).

The volume of all marsh areas at the BBE study site was then calculated using ArcMap 10.3 tools (Supplementary Figure S1), following the approach of Gailis et al. (2021) for BBW (Supplementary Figure S2). The volume was calculated using 140 data points (128 depth profiles and marsh depths obtained from the 12 sediment cores) to interpolate a surface elevation using the Kriging geostatistical method, which uses surface elevation to calculate an estimate of the volume. In the process, the highest surface elevations are made lower than the deepest points recorded (0.98 m) and the lowest surface elevations are then higher than the shallowest points recorded (0.07 m). Uncertainty is introduced into these volume estimates due to the shallow depths recorded (0.07 – 0.98 m) relative to the distance between measurements (25 – 75 m) and the variable topography of the marsh (Amante, 2018; Gailis et al., 2021). Thus, a second, simplified method was used to estimate volume by multiplying the area by the average of the uncompacted core lengths, for the same bounded area estimated via the kriging geostatistical method.

2.2.3 Change in marsh area

Google Earth Pro tools were used with a 50 x 50 m resolution Google Satellite base map (Google Earth Pro, 2019) and an aerial photograph taken on July 28, 1930 (provided by NRCan, National Earth Observation Data) to estimate the extent of the marsh in BBE and MB. The 1930 aerial photograph was overlaid on the Google Satellite base map to quantify changes in the extent of the marsh. This study acknowledges that significant changes to the marsh area were likely made during the construction of Highway 99 in South Surrey between the 1940’s and 1960’s. A similar analysis was previously completed for BBW as part of Gailis et al. (2021), and these results are compared in the Discussion.

Tidal stage and water clarity can be important considerations for delineating marsh extent in aerial photographs and can introduce some uncertainty in delineating marsh edges. The historical aerial photographs we used were generally acquired during lower tidal stages, when the marsh platform and seaward edge were visible. We screened images and selected those with minimal tidal inundation where possible. High marsh vegetation is rarely submerged, while low marsh vegetation may be inundated during high tides but remains partially emergent. Thus, both zones were generally discernible in the imagery, even if water was present. Because we (1) preferably used imagery with lower tidal coverage, (2) cross-checked vegetation signatures across time series, and (3) focused on changes visible at decadal scales, we are confident that the patterns of marsh expansion and retreat we report are robust.

2.3 Field sampling

When defining our high and low marsh sections for C and vegetation analyses, we carefully placed sampling points within the areas of each marsh zone, avoiding the intermediate ecotonal strip where possible. This ensured that the measurements represent the distinct high and low marsh communities. At BBW, 22 sediment cores (13 high marsh, 9 low marsh) were collected between 2014 and 2018 following standard protocols (Gailis et al., 2021). Gailis et al. (2021) measured vegetation composition, marsh volume and area, C stocks and C accumulation rates (from three high marsh and three low marsh cores). For this study, 18 additional cores (8 high marsh, 10 low marsh) were collected in 2020 and 2021 from BBM, BBE, and MB. At BBM, four sediment cores were collected along a transect perpendicular to the dike, with two cores each from the high and low marsh zones. At BBE, three transects of four sediment cores each (12 cores total) were collected perpendicular to the dike with an equal number of cores from the high and low marsh zones on each transect. In MB, a transect was chosen perpendicular to the dike from which two sediment cores were collected, both from low marsh areas (Table 1).

Table 1

Table 1. Sediment core sampling information collected in Boundary Bay, Delta/Surrey, B.C.

Porewater salinities were collected at 22 sampling sites within the four study areas, biweekly from June to August 2020, and monthly from November 2020 to May 2021 (Figure 1). One low marsh porewater collection site was located at MB; two high marsh and two low marsh sites each were located at BBE, BBM, and BBW, with one additional salt panne site at BBW. Porewater was collected using a porewater sipper and a syringe to suction water from the soil at a depth of 25 cm. Approximately 60 mL of soil porewater were extracted, and the salinity was tested using hand-held YSI 556 and YSI ProQuatro conductivity meters (± 0.1 ppt) calibrated for salinity and conductivity. Annual averages along with summer and winter seasonal averages were calculated for each location (Supplementary Table S1).

Cores were collected using an AMS sediment corer pushed into the ground to the depth of refusal. All cores were collected into PVC tubes of 4.7 cm diameter at the field site and stored in the Parks Canada laboratory refrigerator at 4˚C until further analyses. Cores were kept upright from collection time through transport to the laboratory and inside the refrigerator, until they were sectioned. As the cores experienced varying levels of compression during field sampling, a compression factor (CF) was applied to correct for this error, in all cores. The CF was calculated by dividing the depth of core penetration (i.e., length of the hole) by the length of core recovered; this factor was then used to estimate the uncompressed core length.

Sampling included 50 x 50 cm quadrat vegetation surveys at all coring sites (Supplementary Table S2). Additional sample collection in BBE involved the collection of 128 marsh depth profiles (i.e., field measurements of marsh thickness to refusal, used to map marsh depth and volume), and 18 cores collected in 2020 and 2021.

2.4 Laboratory analyses

Sediment cores (n = 18) were sectioned into 1-cm increments in the Parks Canada laboratory (the other 22 sediment cores were already processed, in the same way, for Gailis et al., 2021). Differences in soil composition along the cores were observed in both high and low marsh cores. First, thin layers of white-colored clay mixed with peat or silt were found between 15–62 cm below the surface across four high marsh cores and three low marsh cores from BBM and BBE. The presence of clay material in our marsh cores is consistent with previous findings showing that 15% of the sediment load in the lower Fraser River is clay (Dashtgard and La Croix, 2015). Other sources of clay could be sediment input from local rivers such as Nicomekl and Serpentine, and sediment breaking off from Point Roberts bluffs (Shepperd, 1981; Swinbanks and Murray, 1981). Second, wood debris was common in all four cores at BBM, as well as in two high marsh cores and four low marsh cores in BBE. Visible chunks were removed from the sediment prior to weighing and loss-on-ignition (LOI) measurements. These wood fragments represent allochthonous C, i.e., C derived from outside the marsh ecosystem, in contrast to autochthonous C, which is produced within the marsh by vegetation and microbial activity. Our study focuses on autochthonous C because it reflects the marsh’s intrinsic C sequestration and storage capacity (realizing that marine-derived allochthonous C exists in all salt marsh habitats). Although stable C isotope analysis was not performed in this study, removing large allochthonous materials such as wood helps reduce potential bias in estimating marsh-derived C stocks. Additionally, wood removal was particularly important because this C is already accounted for in terrestrial C budgets. We also removed other obvious non-sediment materials, including shell fragments and animal remains, during processing.

All samples were weighed using a laboratory analytical balance (± 0.01 mg) and then oven-dried at 60˚C for 72 hours. The dry weight (DW, g) was then used to determine the dry bulk density (DBD, g cm^-3) of each sample, relative to the volume ($V_w$, 17.35 cm³) of each sample (Howard et al., 2014).

Samples were ground individually using a 500 mL porcelain lab mortar and pestle. Loss-on-ignition (LOI) was performed on dried sediment samples in two steps to estimate organic and inorganic carbon content. First, ground soils were combusted in a muffle furnace (4 hrs at 550 ˚C) to burn off all the organic compounds (Howard et al., 2014). The samples were then weighed to calculate the LOI (%) (Equation 1):

$\begin{array}{ll}\% LO{I}_{550}=\left(\frac{D{W}_i-D{W}_f}{D{W}_i}\right)\times 100& (1)\end{array}$

where $D{W}_i$ is the initial dry weight and $D{W}_f$ is the dry weight after burning.

Second, samples were heated at 1000˚C to remove carbonates, allowing distinction between organic and inorganic C fractions. To quantify the fraction of organic C, a subset of samples (n = 39) was burned again in the muffle furnace for 2 hours at 1000˚C to determine its inorganic C (IC) content (Heiri et al., 2001) (Equation 2):

$\begin{array}{ll}\%LO{I}_{1000}= \left(\frac{D{W}_{550}- D{W}_{1000}}{D{W}_i}\right) \times 100& (2)\end{array}$

where $D{W}_i$ is the initial dry weight, $D{W}_{550}$ is the dry weight after the 550˚C burn and $D{W}_{1000}$ is the dry weight after the 1000˚C burn.

The percentage of inorganic C (IC) was assumed to be negligible in all samples, based on previous coulometric measurements from western Boundary Bay, which showed that IC content to be less than 0.01% (Gailis et al., 2021). Due to inaccessibility of laboratories during the COVID-19 pandemic, an elemental analysis could not be completed to determine an in situ empirical relationship between %C and %LOI for the eastern portions of Boundary Bay. Given the geographical proximity of our sites, we applied the regression equation established by Gailis et al. (2021) for western Boundary Bay (r² = 0.97) to calculate %C from %LOI for all mid, eastern and Mud Bay sites (Equation 3).

$\begin{array}{ll}\%C=0.44(\%LOI)-1.3& (3)\end{array}$

The regression produced a small negative intercept, which can result in negative %C values for low-LOI sediments. Several previous studies have indicated that LOI has the potential to overestimate soil C because the ignition process drives off both organic matter as well as water bound in any clay minerals that are present in the sample (e.g. Howard, 1966; Howard and Howard, 1990; Santisteban et al., 2004). To preserve the empirical relationship observed in the calibration dataset, the intercept was not forced through zero. Any negative %C values were subsequently capped at zero, as negative C is physically impossible. This approach allows realistic C stock estimates while maintaining the integrity of the LOI–%C relationship.

2.5 Carbon stocks and C storage

2.5.1 Soil carbon densities and carbon stocks

Carbon stocks were quantified by measuring the soil C density (SCD) for each 1-cm sample down to the basal peat layer (n = 18, eight for high marsh, 10 for low marsh, measured to DoPs: depth of peat layer). We measured C stocks down to the base of peat layer due to the shallow depth and young age of BC salt marshes (Chastain et al., 2022). For comparison, average marsh depths were 23.4 cm in nearby Clayoquot Sound along the southwest coast of Vancouver Island, approximately 49.17°N, 125.93°W and 23 cm in western Boundary Bay. Basal peat ages were between 13 and 140 yrs (Clayoquot Sound) and between 37 and 67 yrs (western Boundary Bay) (Gailis et al., 2021; Chastain et al., 2022).

The base of the peat layer was not reached in two cores: BBM1H1 and BBE2H2. We still include these cores in our calculations of uncompressed peat layer depths and C stocks and note here that our inability to reach the base of the peat layer results in an underestimation of the peat layer depth and C stocks for the BBM and BBE high marsh zones.

Following the approach of Chastain et al. (2022), we used the age estimates to calculate the 30-year C stocks in a subset of cores with ²¹⁰Pb dating. We chose the 30-year C stock because (a) this time horizon allows us to compare C stocks across regional studies (e.g. Chastain et al., 2022); and (b) a 30-year time horizon is reasonable both for policy considerations and for comparing with the C sequestration impacts of any recent restoration projects that have taken place within the last 30 years.

Soil C density (SCD, ${\text{g C cm}}^{-3}$) was derived from the measured dry bulk density (DBD) and percentage of C (%) for each centimeter interval sampled (Gailis et al., 2021; Howard et al., 2014) (Equation 4):

$\begin{array}{ll}SCD\left(\frac{g C}{c{m}^3}\right)= \left(\frac{\%C}{100}\right) \times DBD \left(\frac{g}{c{m}^3}\right)& (4)\end{array}$

The total C stock for each core (g C cm^-2) was calculated by summing all SCDs at 1-cm intervals down to the base of peat layer in each core (Chastain et al., 2022) (Equation 5).

$\begin{array}{ll}Cstoc{k}_{core} \left(\frac{g C}{c{m}^2}\right)= {\sum }_{i= 0 }^{n}SC{D}_i \times 1 cm& (5)\end{array}$

where $n$ is the number of slices within the core (cm) and $SC{D}_i$ is the SCD of each 1-cm interval (g C cm^-3).

2.5.2 Compression factors

As the cores experienced various levels of compression during field sampling, a compression factor was used to compensate for the error. To calculate the compression factor (CF) (Table 1), the length of core penetration was divided by the length of core sample recovered, and this was used to estimate an uncompressed core length (Equations 6, 7):

$\begin{array}{ll}CF= \frac{length of core penetration (cm)}{length of core (cm)}& (6)\end{array}$

$\begin{array}{ll}Uncompressed Length (cm)=CF x Compressed Length of Core (cm)& (7)\end{array}$

To calculate the uncompressed SCD, the compressed SCD for each 1-cm soil section was divided by the compression factor for that specific core (Howard et al., 2014) (Equation 8):

$\begin{array}{ll}Uncompressed SCD \left(\frac{g C}{c{m}^3}\right)=\frac{ Compressed SCD \left(\frac{g C}{c{m}^3}\right)}{CF} & (8)\end{array}$

The compression factor is applied to the entire length of the core, assuming that all sections are compacted equally. However, this assumption is likely an oversimplification, as factors such as soil type, moisture content, and dry bulk density – which influence compression – are likely to vary throughout the core (Howard et al., 2014; Morton and White, 1997).

2.5.3 Marsh carbon storage

We calculated total C storage (Mg C) for high and low marsh zones separately to avoid generalizations between high and low marsh C stocks, SCDs, and depth of the basal peat layers (DoPs). We combined the total C storage for high and low marsh zones by adding their values to obtain the total marsh C storage down to the base of the peat layer. Then, we calculated C storages at each of the four study sites separately and add them together to derive the total C storage of Boundary Bay marsh. Although previous C storage was calculated for BBW by Gailis et al. (2021), we have re-calculated C stocks to the basal peat layer for each BBW core to maintain consistency with calculations from the eastern portion of the bay. Based on this recalculation, we also revise the BBW C storage, as our areal estimates for the BBW marsh are also larger than reported in Gailis et al. (2021).

Three methods were used to calculate total C storage. The traditional method estimates total C storage by calculating the average core C stocks (Mg C ha^-1) to the base of the peat layer in each core, using compressed SCDs and multiplying it by the area (ha) of the marsh (Gailis et al., 2021). We used compressed SCDs in this calculation because the amount of C down to the base of peat layer of each core remains the same regardless of whether the core is compressed (Equation 9), as long as the compressed depth (to peat layer) is used with the compressed SCD.

$\begin{array}{ll}Cstorag{e}_{Marsh }(Mg C )= A\times {\overline{Cstock}}_{core} & (9)\end{array}$

where ${\overline{Cstock}}_{core}$ is the average core C stock to the base of the peat layer calculated using compressed SCDs, and $A$is the total area of the marsh.

The other two methods, kriging volume and simplified volume methods, use the volume of the marsh (m³) to calculate the total C stock.

The simplified volumetric approach multiplies the marsh area by mean uncompressed core length (averaged from 12 cores for BBW, 4 cores for BBM, 12 cores for BBE, 2 cores for MB, total n = 30) down to the peat layer to derive the simplified volume (m³), which is then multiplied by the average uncompressed SCD (estimated for each of the four study sites separately). This calculation is done for each of the four study sites separately and then added together to derive the Boundary Bay marsh total (Equation 10):

$\begin{array}{ll}Cstorag{e}_{Marsh }(Mg C)= A\left(m^2\right) \times Z_{core} \times Average Uncompressed SCD& (10)\end{array}$

where $A$ is the marsh area for each of the four study sites and $Z_{core}$ is the average uncompressed core length (m) down to the peat layer within each of the four study sites.

For the third approach, the kriging volume of the marsh (m³) is calculated using the kriging geostatistical tool on ArcMap 10.3 derived from depth profiles. We note that, using this approach, the marsh volume is calculated using the total thickness of marsh deposits, from the surface to refusal (measured using depth profiles). Thus, this measurement includes peat and any underlying mineral layers to refusal. This method was only completed for our two sites with extensive depth profiling: BBW (n = 176; Gailis et al., 2021) and BBE (n = 140). Although this measurement does not estimate C storage for the entire Boundary Bay marsh, it is a useful reference for comparing C storage to other estimation methods for two of our most extensively cored study sites. The kriging volume for the BBW and BBE is then multiplied by the average of uncompressed SCDs for both sites to derive the marsh C stock (Gailis et al., 2021) (Equation 11):

$\begin{array}{ll}Cstorag{e}_{Marsh }(Mg C )= V_{kriging}\times Average Uncompressed SCD& (11)\end{array}$

where $V_{kriging}$ is the geostatistical kriging-derived marsh volume (m³).

2.6 Carbon accumulation rates

Carbon accumulation rates (CARs) were estimated based on the distribution of the natural radionucide ²¹⁰Pb (Appleby, 2001). Subsamples from four BBE cores (BBE1H1 and BBE3H1 for high marsh, BBE1L and BBE3M for low marsh) were sent to Flett Research Ltd (Winnipeg, MB) for ²¹⁰Pb radiometric dating analysis. Cores were selected to best represent the spatial variability and typical sediment layers in the marsh. All cores were corrected for depth compression using the compression factors calculated above to ensure accurate sedimentation rates (Morton and White, 1997; Gailis et al., 2021).

All samples sent for radiometric dating were freeze-dried for four days, and the dried material was ball-milled prior to radioisotope analyses. Between 7 to 13 dried subsamples per core were analyzed for ²¹⁰Pb content. Activities of ²¹⁰Pb (Bq kg^-1) were determined by α-spectrometry through its granddaughter product ²¹⁰Po, assumed in secular equilibrium. The atmospheric, unsupported, or excess ²¹⁰Pb (²¹⁰Pb_xs) fraction was used to derive the age-depth model and was determined as the difference between the total ²¹⁰Pb activity and its parent nuclide ²²⁶Ra activity (Figure 2; Supplementary Table S4) (Eakins and Morrison, 1978; Mathieu et al., 1988). One to three ²²⁶Ra measurements per core were made and assumed to equal the supported ²¹⁰Pb (²¹⁰Pb_sup) (Arias-Ortiz et al., 2018; Callaway et al., 2012; Appleby and Oldfield, 1978; Ritchie and McHenry, 1990). The ²²⁶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 days. Samples were sealed and stored for at least three weeks before counting to ensure secular equilibrium of ²²⁶Ra daughters. The ²²⁶Ra content was determined by counting ²²²Rn activity for 60,000 seconds (Minimum Detectable Activity (MDA) = 0.0167 Bq kg^-1) (Durham and Joshi, 1980; Mathieu et al., 1988; Pennington et al., 1973; Anderson et al., 1987).

Figure 2

Figure 2. Natural log-transformed excess ²¹⁰Pb activity (ln(210Pb_xs)) plotted against cumulative mass for high marsh cores (a–e) and low marsh cores (f–j) from all Boundary Bay sites. Dashed lines show the linear regressions.

As documented previously in Gailis et al. (2021), ²¹⁰Pb and ²²⁶Ra analyses were also conducted on six sediment cores from western Boundary Bay (Table 2). Four cores were sent to GEOTOP laboratories at the Université du Québec à Montréal (Montreal, QC), and two cores were sent to Flett Research Ltd (Winnipeg, MB). As with the cores in eastern Boundary Bay, all samples sent for radiometric dating were oven-dried (72 hrs at 60 °C), weighed to determine DBD, and homogenized using a mortar and pestle. Between eight and 11 dried sub-samples per core were analyzed for ²¹⁰Pb measurements. One ²²⁶Ra measurement per core was also conducted.

Table 2

Table 2. Geochronological, carbon, and sediment depth information used to estimate individual core as well as high and low marsh means ± SD of sedimentation rates (SAR), mass accumulation rates (MAR), carbon accumulation rates (CAR), and 30-year C stocks across all of Boundary Bay.

Sediment chronologies for all cores were estimated using the Constant Flux: Constant Sedimentation (CF: CS) model (Krishnaswamy et al., 1971). The CF: CS model assumes that the supply of ²¹⁰Pb_xs to the sediment surface and the mass accumulation rate (MAR, g m^-2 yr^-1) remain constant through time. Under these assumptions, the MAR can be obtained from the slope (cm² g^-1) of the linear best-fit line of the relationship between the natural log of the ²¹⁰Pb_xs activity (Bq kg^-1) against the cumulative mass (g cm^-2) using Equation 12:

$\begin{array}{ll}\text{MAR} (g c{m}^{-2 }y{r}^{-1})=\frac{0.0311 y{r}^{-1}}{slope (c{m}^2 g^{-1})}& (12)\end{array}$

where 0.0311 yr^-1 is the radioactive decay constant of ²¹⁰Pb.

Because changes in accumulation rate may have occurred, the CF: CS model was applied in a piecewise way in those core intervals with differing sedimentation rates (Figure 2).

Carbon accumulation rates (CARs) were calculated for the whole core, by multiplying the average uncompressed percent carbon (%C_uncompressed) to the base of the peat layer by the MAR (Equation 13):

$\begin{array}{ll}CAR (\text{g C }{m}^{-2}y{r}^{-1})=\frac{\text{\%C}}{100}\times MAR (\text{g }{m}^{-2}y{r}^{-1})& (13)\end{array}$

2.7 Statistical analysis

All data were tested for normality using the Shapiro-Wilk test for normality. T-tests were conducted for DBDs, SCDs, %C, C stocks, C storages, MARs, and CARs to test for any significant differences between high and low marsh sites, and between the four study areas (BBW, BBM, BBE, MB). The significance level of all the tests was set at p = 0.05. All statistical analyses were performed in R (RStudio, 2015).

All raw data for BBW were taken from a previous study by Gailis et al. (2021). DBDs, SCDs, %C, C stocks and C storages were recalculated in this study down to the base of the peat layer to ensure consistent treatment of data at all sites. The t-tests and plots provided in this study are based on the recalculated values of DBDs, SCDs, %C, C stocks and C storage in BBW.

3 Results

3.1 Sediment properties

Uncompressed full lengths of the collected marsh cores ranged from 23 to 78 cm for all of Boundary Bay (Table 3). Compression occurred in all cores during field sampling and was similar in both high and low marsh. Compression factors ranged between 1.00 and 3.00 in all cores (Table 1, n = 39, not including W-A-H sediment pit, which was dug for Gailis et al. (2021) and experienced no compression). W-A-H was a sediment pit that was dug until sand was reached, to better understand the degree of compression and to explore what sedimentary layers occurred below the depth of refusal.

Table 3

Table 3. Summary of uncompressed core sediment data (± SD) for depth of core, depth of peat layer, dry bulk density (DBD), average soil carbon density (SCD), average %C, and C_stock down to the base of peat layer for cores (with documented compression factors) in Boundary Bay, Delta/Surrey, B.C.

Our porewater salinity measurements indicate that the marsh is largely polyhaline (e.g. sits above a critical salinity threshold of 18) in BBW and BBM, but mesohaline in BBE and MB (salinities between 5 and 18). Complete year-round porewater salinity measurements were obtained at several sites across Boundary Bay (Supplementary Table S1). At BBW, both high (BB2-5) and low marsh (BB2-4) sites consistently measured salinities above 18 ppt in both summer and winter months, 100% of the time. Moving eastward to BBM, low marsh sites BBM1M and BBM1L measured salinities above 18 ppt 80% and 60% of the time, respectively. Further east, high marsh sites BBE1H2 and BBE1H1 measured salinities above 18 ppt only 20% of the time.

In the eastern portion of Boundary Bay, soil horizons in high and low marsh cores differed in their composition (Supplementary Figure S3). High marsh cores generally had a thick top organic peat layer composed of dark brown organic material (uncompressed average ± standard deviation (SD) = 29 ± 28 cm), followed by layers with a mix of peat and silt or silty sand. Cores BBE2H2 and BBE2H1 had the thickest (uncompressed) organic peat layers at 66 and 63 cm, respectively. The bottom of two high marsh cores consisted of sand material, two consisted of clay, two ended in silt, one in peat, and one in peat followed by wood. In contrast, low marsh sediment cores had thinner top organic layers compared to high marsh cores (uncompressed average ± standard deviation (SD) = 18 ± 13 cm). The bottoms of all low marsh cores ended in a sand layer.

In all cores, the percent weight of C (%C) was highest in the top peat layers and declined at deeper depths (Figure 3). Maximum (uncompressed) values of %C approached 25% in the peat layers of high marsh cores (BBE3H2 from BBE, BB1 from BBW) and exceeded it in one of the low marsh cores (BBM1L from BBM). Bottom sand layers contained the lowest %C, reaching 0% C in both high and low marsh cores.

Figure 3

Figure 3. Changes in uncompressed percentage of carbon (%C) vs uncompressed sediment depth (cm) for (a) BBW high marsh, (b) BBW low marsh, (c) BBM high marsh, (d) BBM low marsh, (e) BBE high marsh, (f) BBE low marsh, and (g) MB cores. BBW data are from Gailis et al. (2021); all other data were collected for this study.

Across all Boundary Bay marsh cores (Table 3; Figures 3a–g, 4a–f), the average uncompressed %C ± SD down to the base of the peat layer was significantly higher in the high marsh (9.6 ± 4.1%; n = 16) than in low marsh cores (5.4 ± 3.3%; n = 15) (p<0.005). Similar patterns were observed in BBE (high marsh = 9.8 ± 5.9% (n=6), low marsh = 3.8 ± 1.7% (n = 6)) (p<0.05). In BBW, differences were observed, but did not meet the significance threshold (high marsh = 9.4 ± 3.3% (n = 8), low marsh = 6.2 ± 2.4% (n = 5), p = 0.08). No significant difference in %C was observed between high and low marsh in BBM (high marsh = 9.4 ± 1.5% (n = 2), low marsh = 10.8 ± 1.1% (n = 2) (p>0.05). MB had the lowest overall %C values at 2.7 ± 1.3% (n = 2) (Figure 3g), although slightly higher average %C values were observed in the low marsh cores of BBE (3.8 ± 1.7%) and BBW (6.2 ± 2.4%).

Figure 4

Figure 4. Box and whisker plots showing distribution of average (uncompressed) percent carbon (%C) in each core down to the base of peat layer for each marsh region. The middle line is the median and the top and bottom of the box are quantiles (Q1 and Q3), and the error bar is the largest and smallest value (dots are outliers). Data are shown for high and low marsh zones in (a) BBM (b) BBE, (c) MB, (d) BBW, and (e) all Boundary Bay cores. All sites (high and low marsh data combined) are compared in (f). T-tests show significant differences (indicated with asterisk) between low and high marsh %C for (b) BBE (t-value=2.6033, p-value=0.0263, p<0.05), and (e) all Boundary Bay combined (t-value=3.3670, p-value=0.0022, p<0.01), but not for (a) BBM (t-value=1.0985, p-value=0.3866, p>0.05) or (d) BBW (t-value=1.9051, p-value=0.0832, p<0.1).

As with the %C content, the uncompressed SCD decreased from the surface down to the base of the peat layer (Supplementary Figure S5), and the uncompressed SCDs were also significantly higher in the high marsh (0.028 ± 0.008 g C cm^-3; n = 15) than in the low marsh zones (0.020 ± 0.008 g C cm^-3, n = 13) across all study sites in Boundary Bay (p<0.05) (Supplementary Figure S6). No significant differences in SCDs were observed between high and low marsh cores in BBE, BBM, and BBW (p>0.05). Uncompressed SCDs in the western portion of the marsh (BBW, BBM, n = 15) were significantly higher than that of eastern study sites (BBE, MB, n =13) at Boundary Bay (p<0.05).

The uncompressed DoP was recorded to approximate the base of the marsh (Table 3). Across Boundary Bay, no significant differences were observed in DoP between low (19 ± 13 cm) and high (29 ± 20 cm) marsh areas (p>0.05) (Table 3; Supplementary Figure S6). The maximum uncompressed DoP was recorded in core BBE2H2 at 66 cm, and the minimum was recorded in cores BBE1L and MB1L at 2 cm each. No significant difference was observed in DoP between cores from the western and eastern portions of Boundary Bay (p>0.05).

3.2 Carbon stocks, marsh area, carbon storage, and accumulation rates

3.2.1 C stocks

Overall, C stocks down to the base of the peat layer for Boundary Bay marsh averaged 54 ± 36 Mg C ha^-1 (n = 39; Table 3). Among all sites, BBM (n = 4) had significantly higher C stocks than all other sites at Boundary Bay (n = 28) (Figures 5a–f) (p<0.01). In addition, C stocks in the western sites (BBW, BBM, n = 18) were significantly higher than in the eastern sites (BBE, MB, n = 14) (p<0.01) (Figure 5f).

Figure 5

Figure 5. Distributions of average core C stocks (Mg C ha^-1) down to the base of peat layer for high and low marsh cores for marsh sites and all of Boundary Bay combined. Areas sampled include (a) BBM (b) BBE (c) MB, (d) BBW, (e) All Boundary Bay, and (f) all C stock estimates (both high and low marsh) separated by site. Differences between high and low marsh C stocks are not statistically significant for any marsh site but are statistically significant (t = 2.2482; p < 0.05) when all sites are combined as high vs low (panel e, indicated with asterisk). The middle line is the median and the top and bottom of the box are quantiles (Q1 and Q3), and the error bar is the largest and smallest value. Dots represent outliers.

When considering all sites, high marsh C stocks at Boundary Bay (71 ± 37 Mg C ha^-1) were significantly higher than low marsh C stocks (41 ± 36 Mg C ha^-1) (Figure 5e) (p<0.05). In BBE, average C stocks ranged from 21 to 112 Mg C ha^-1 (average = 55 ± 41 Mg C ha^-1) in high marsh and 0 to 40 Mg C ha^-1 (mean = 29 ± 16 Mg C ha^-1) in low marsh (Figure 5b, Table 3). The average of 27 ± 15 Mg C ha^-1 for low marsh cores is skewed by the negligible C stock of core BBE1L (0 Mg C ha^-1). When BBE1L is not considered, the average low marsh C stock is substantially higher at 34.5 ± 7.4 Mg C ha^-1. MB reported the lowest average C stock for low marsh cores at 6.5 ± 8.8 Mg C ha^-1 (n = 2) (Figure 5c).

3.2.2 Marsh area

The total area of the Boundary Bay salt marsh is 222 ha, with high and low marsh each accounting for 50% of the total marsh area. In the western portion of the Boundary Bay marsh, BBW covers a total area of 80 ha, with high marsh accounting for 62.5% (50 ha) of the total marsh area. BBM contains the greatest marsh area among all sites at 84 ha, with high marsh accounting for 46% (39 ha) of the total marsh area. In the eastern portion of the Boundary Bay marsh, BBE spans 52 ha, with high marsh accounting for 42% (22 ha) of the total marsh area. MB marsh spans an area of 6 ha, with low marsh accounting for 100% of the total marsh area.

Comparisons of our marsh area estimates with aerial photographs from 1930 indicate that both the eastern and western portions of the marsh have changed in area over the past ~90 years (Figures 6a, b). Previous work by Gailis et al. (2021) indicates that the western portion of the Boundary Bay marsh (BBW and BBM) expanded by about 26 ha (~20%) between 1930 and 2018. In contrast, Figures 6a, b show that the eastern portion of the Boundary Bay marsh (BBE, MB) was reduced by about 18 ha (~35%) in the past 91 years (1930 - 2021) (Figure 6b).

Figure 6

Figure 6. Comparison of 1930 and 2021 marsh extent for BBE and MB in eastern Boundary Bay, BC. (a) Air photo of BBE and MB from 1930 was superimposed onto Google Earth 2021 base map. Yellow line represents present-day dike position. Black arrow indicates a larger than present-day channel in BBE high marsh area. Red circles indicate marsh areas cleared out and filled in during Highway 99 and dike construction post-1930s. Base Map Source: NRCan, National Earth Observation Data. (b) Map of present-day BBE and MB relative to the leading edge of marsh in 1930. White line represents the marsh edge in 1930 obtained from panel (a), and yellow line represents present-day dike position. Base Map Source: Google Earth Pro 2021. In both panels, red dots represent the four dated cores from eastern Boundary Bay and white dots represent all other sediment cores.

Several changes have occurred in the eastern portion of Boundary Bay since 1930. First, a large channel (indicated by the black arrow in Figure 6a) appears to have inundated a larger portion of the BBE high marsh area in 1930. Additionally, the construction of Highway-99 in South Surrey in 1942 appears to have cleared and filled some areas of the low marsh located between BBE and MB sites (indicated by red circles in Figure 6a), likely resulting in a smaller, present-day marsh. Dike renovations also appear to have impacted the shape and size of the MB marsh, primarily near the coring sites. The MB marsh area on which cores were extracted appears to have developed after 1930 (Figure 6b).

3.2.3 Marsh C storage

Total C storage was estimated using three methods for sites BBW (Gailis et al., 2021) and BBE: the traditional method, the simplified volumetric method, and the kriging volumetric method (Supplementary Table S5). At these locations, the total C storage estimates made using all three approaches were not statistically different. At BBW, the kriging volumetric method estimated total C storage of 5,880 ± 1,680 Mg C, which is within the uncertainty of estimates obtained using both the simplified volumetric approach (5,320 ± 1,520 Mg C) and the traditional method (5,150 ± 1,470 Mg C). For BBE, the three approaches also produced estimates that were not significantly different, given their uncertainties (2,520 ± 960 Mg C for kriging volume; 2,510 ± 960 Mg C for simplified volume; and 2,620 ± 1,000 Mg C for the traditional method). We have chosen to report the estimates made from the simplified volumetric approach (although all calculations are provided in Supplementary Table S5).

The total estimated C storage for all Boundary Bay marsh was 17,360 ± 4,960 Mg C. Of this C, 61% (10,540 ± 2,720 Mg C) was found in the high marsh zone. The remaining 39% of C was found in the low marsh zones (Supplementary Table S5). The combined C storage in western study sites (BBW, BBM) makes up nearly 85% of the total C storage at Boundary Bay marsh. These estimates (and all estimates in this section) reflect the amount of C stored down to the base of the peat layer.

When we consider site-specific differences in C storage, we find that BBM had significantly higher C storage compared to all other sites (p<0.01). Total C storage at BBM was 9,300 ± 1,240 Mg C. Here, the high marsh zone made up only 45% of the total C storage, with estimates of 4,200 ± 560 Mg C. The next highest C storage was found in BBW (5,320 ± 1,520 Mg C). BBW holds the largest amount of C in the high marsh zone (82%). (Supplementary Table S5). BBE had a total C storage of 2,520 ± 960 Mg C. Here, 62% of the total C storage was found in the high marsh zone (Supplementary Table S5). The least amount of C was estimated for MB, which was two orders of magnitude lower, with a total C storage of 26 ± 29 Mg C. This site had no high-marsh sediments, and so all C was found in the low-marsh zone.

3.2.4 C accumulation rates and total marsh annual C sequestration

For all of Boundary Bay, high marsh MARs ranged from 431 ± 22 to 1,634 ± 134 g m⁻²y⁻¹ (average = 922 ± 506 g m⁻²y⁻¹), and low marsh MAR ranged from 757 ± 69 to 1,745 ± 235 g m⁻² y⁻¹ (average = 1,255 ± 464 g m⁻²y⁻¹) (Table 2). No significant differences were found between the high and low marsh MARs (p>0.05).

High marsh CARs ranged from 24 to 152 g C m⁻² yr⁻¹ (average = 89 ± 56 g C m⁻² y⁻¹). Low marsh CARs ranged from 36 to 105 g C m⁻² yr⁻¹ (average = 69 ± 32 g C m⁻² y⁻¹). No significant differences were found between the high and low marsh CARs (p>0.05) (Table 2). CARs averaged 80 ± 45 g C m⁻² yr⁻¹ for all sites at Boundary Bay. When extrapolated to the entire 222 ha of marsh, this equates to an annual C sequestration rate of 178 ± 100 Mg C yr^-1.

The basal ages of the sediment cores (depth of basal peat layer) ranged from 26 to 89 years for all cores in Boundary Bay (Table 2), with no statistically significant difference between the high and low marsh settings sampled (p>0.05).

The CF: CS model could not be applied to the low marsh core BBE1L, as it did not meet the following assumptions: sufficient core length to capture the full inventory of atmospherically derived ²¹⁰Pb, and a constant sediment accumulation rate. The irregular shape of BBE1L ²¹⁰Pb profile indicates intense mixing, erosion, and highly variable sediment accumulation rates at this sampling site.

The 30-year C stocks were calculated for the nine cores where ²¹⁰Pb dating was successfully applied. These stocks ranged from 12 to 69 Mg C ha^-1, with an average of 32 ± 16 Mg C ha^-1 across all of Boundary Bay. No statistically significant differences were found between the five high marsh cores (average = 36 ± 20 Mg C ha^-1) and four low marsh cores (average = 27 ± 11 Mg C ha^-1) (p>0.05).

4 Discussion

4.1 Carbon stocks

Higher C stocks are buried in the western portion of Boundary Bay (BBM, BBW) compared to the east (BBE, MB), in part due to high C stocks in both the high and low marsh zones in BBM (Figure 4f). Compared to eastern sites such as BBE and MB, we found that western Boundary Bay sites, especially BBM, have a more mature and diverse plant canopy that enables autochthonous C build-up, along with allochthonous organic matter brought in by tides to be trapped more easily in the soil. Our field observations indicated that substantial amounts of eelgrass rack were deposited in western Boundary Bay over the winter months. This above-ground accumulation likely leads to thicker organic matter layers at the top of the sediment cores, resulting in higher C stocks compared to other parts of the marsh (e.g. Granse et al., 2024).

Unlike all other parts of Boundary Bay, no significant difference in C stocks was observed between high (108 ± 11 Mg C ha^-1) and low marsh zones (115.0 ± 4.3 Mg C ha^-1) at BBM. The relatively high C stocks in the BBM low marsh could be attributed to the somewhat bay-like shape morphology of this location (Figure 1), which facilitates the accumulation of higher-than-normal amounts of seagrass and other organic matter (Dashtgard, 2011). We acknowledge, however, that fewer cores were collected in BBM (n = 4 covering 84 ha) compared to BBW (n = 22 covering 80 ha). Additional cores along multiple transects within BBM might reveal greater variability in C stocks, similar to patterns observed in other areas of the Boundary Bay marsh.

Lower C stocks in the eastern portions of Boundary Bay (BBE, MB) compared to the western areas (BBW, BBM) could be attributed to differences in sedimentation patterns influenced by coastal geography, sediment sources, and circulation dynamics (Swinbanks and Murray, 1981). Our study sites extend across the 19 km length of Boundary Bay, and each area of the marsh has differences in geographical placement, tidal inundation, vegetation, sedimentation levels and C accumulation. The western portion of the marsh is likely expanding due to the accumulation of sediment breaking off the Pleistocene cliffs at Point Roberts (Swinbanks and Murray, 1981). The abundance of nooks and sheltered “bay-like” areas in mid and western marsh areas also facilitates the accumulation of sediments and organic debris. Simultaneously, the eastern portions are likely receding due to the river and water currents, as well as wind dynamics in the Bay (Dashtgard, 2011; Swinbanks and Murray, 1981). While western Boundary Bay is somewhat protected from wind and ocean currents by the Point Roberts peninsula, eastern Boundary Bay is more directly exposed to high waves and strong ocean currents that could lead to mechanical erosion and the hummocky marsh edges that are observed in the eastern part of the marsh. These differences lead to the observed variability in C stocks across the marsh.

Mud Bay consistently had %C, SCD and C stock values that were lower than the other parts of the marsh, likely due to its thin organic layer (< 6 cm). According to the 1930 satellite photo from NRCan (Figure 6a), all cores at Mud Bay were extracted on a new marsh area formed after 1930 (< 91 years old), likely on areas cleared and filled in during the construction of Highway-99 adjacent to Mud Bay in 1942 (Figure 6b) and dike renovations in 1948. The overall reduction in marsh area in Mud Bay in the past 91 years (1930 - 2021) could be attributed to disturbance through highway construction. We recognize, however, that only two cores were collected from Mud Bay, and that core locations to the east or west of the current coring sites could be older (> 91 years old) and therefore may hold more carbon. Collecting additional samples could result in higher estimates of C stocks at MB if, for example, other places in MB have higher C densities and C stocks than in the two collected cores (the low marsh C stocks in BBE are 4 times higher than in MB). However, the overall small area of MB (5.9 ha of 222 ha) suggests that its contribution to the total C storage of Boundary Bay marsh would remain small: even if the C storage of MB were 4 times higher, the total C stock of MB would still make up less than 1% of all C stored in Boundary Bay.

C stocks in Boundary Bay (41 ± 36 Mg C ha^-1 for low marsh and 71 ± 37 Mg C ha^-1 for high marsh) represent only 16-28% of the global average estimate of 250 Mg C ha^-1 (Chmura et al., 2003; Pendleton et al., 2012). However, our estimates are quite similar to more recent estimates made by Maxwell et al. (2023), who calculated globally averaged C stocks of 79.2 ± 38.1 Mg C ha⁻¹ in the top 30 cm of global salt marshes. The similarity with the Maxwell et al. (2023) estimate makes more sense given that the average peat depth in Boundary Bay, corrected for compression, was only 19 ± 13 cm. We also note that the average uncompressed SCD across Boundary Bay (0.024 ± 0.006 g C cm^-3) is 38% lower than the global estimate of 0.039 ± 0.003 g C cm^-3 (Chmura et al., 2003). This difference may reflect lower C contents in Boundary Bay sediments but also may be influenced by methodological differences in estimating %C from %LOI. In this study, we used a regionally derived empirical relationship between %C and %LOI to calculate the fraction of C in Boundary Bay sediments (Gailis et al., 2021). In contrast, Chmura et al. (2003) used a generalized polynomial relationship from Craft et al. (1991). Notably, Chastain et al. (2022) noted that the use of the Craft et al. (1991) equation resulted in C stock estimates approximately 30% higher than those based on locally derived relationship. Thus, although methodological differences contribute to some uncertainty in C stock estimates, we anticipate that the shallow depth of Boundary Bay marsh is the primary factor.

Comparison with compilations from other parts of the Pacific Coast of North America also suggest that the low C stocks observed in our study region (and other parts of British Columbia) may be unique. In the Pacific Northwest region of the United States (PNW), Kauffman et al. (2020) quantified average total ecosystem C stocks of 417 ± 70 Mg C ha^-1 and 551 ± 47 Mg C ha^-1 for low and high marsh, respectively, but these calculations included C down to depths of 3 m. The top 1 m of soil only accounted for 48-53% of the total C stock. Therefore, global estimates down to 1 m of soil may greatly underestimate C stocks in other parts of the Pacific Coast (Kauffman et al., 2020). We note, however, that C stocks in the surface 30 cm of PNW marshes (88 ± 10 Mg C ha^-1) are comparable to the average C stocks for high marsh sediments in Boundary Bay (71.4 ± 37.4 Mg C ha^-1), where the base of the peat layer averages 29 ± 20 cm. Once again, this comparison suggests that the depth of the Boundary Bay marsh is the cause of its lower C stocks.

When averaged over both high and low marsh areas, peat layer depths at the Boundary Bay salt marsh are relatively shallow, with an average of 24 ± 17 cm. These shallow depths may simply reflect a younger marsh system, with basal peat ages ranging from 26 to 89 years old. The young marsh age can be linked to the relatively recent formation of the Fraser River Delta over the past 5000 years through fluvial sediment deposition. Prior to the formation of the bay-like structure bound by Point Roberts to the west, the exposed Boundary Bay region was more heavily influenced by wind-driven ocean currents compared to present day, which likely hindered marsh establishment and growth (Dashtgard, 2011). In addition, the eastern portion of the bay has been subject to persistent erosion throughout the past 4000 years and lacks sufficient sediment supply to support vertical accretion (Shepperd, 1981). Sedimentological studies also suggest that the current marsh is underlain by a sandy gravel beach-ridge complex (Engels and Roberts, 2005; Dashtgard, 2011). As a result, extrapolating salt marsh C stocks to a uniform 1 m – as is commonly done in global estimates - can substantially overestimate regional C storage, particularly in systems like Boundary Bay where the onset of peat accumulation occurs well above the 1 m reference horizon (Gailis et al., 2021; Chastain et al., 2022).

At the same time, extrapolating to a reference horizon of 1 meter provides a useful estimate of the C storage potential of the Boundary Bay marsh. A simple back-of-the-envelope calculation, based on the average uncompressed SCD (0.025 g C cm^-3) and the area of Boundary Bay marsh (222 ha), suggests that if the marsh were to accumulate 1 meter of sediment, it could store approximately 55,500 Mg C, more than 3 times the current estimated C storage of 17,360 Mg C. This type of calculation holds substantial uncertainties, as it assumes that the marsh accumulation would not be hindered by factors such as storm events (tsunamis), anthropogenic changes, and sea level rise. However, the future C storage potential of Boundary Bay could substantially exceed its present-day C content.

Another approach to understanding C sequestration in Boundary Bay is to consider the 30-year C stocks in cores where ²¹⁰Pb dating has been applied. The 30-year C stocks ranged from 14 to 69 Mg C ha^-1 and averaged 34 ± 16 Mg C ha^-1 across all of Boundary Bay. This value appears smaller (although statistically not different) from the average 30-year C stock calculated for Clayoquot Sound (54 ± 5 Mg C ha⁻¹). Given the 222 ha of Boundary Bay marsh, we can estimate that ~7550 ± 3550 Mg C has accumulated within Boundary Bay over the past 30 years, or 252 ± 118 Mg C yr⁻¹. For comparison, this annual C uptake is roughly equivalent to the annual emissions from 202 passenger vehicles, assuming the average vehicle gets 22.2 miles per gallon (10.6 liters per 100 km), travels 11,500 miles (18,510 km) per year, and every gallon of gasoline burns 8,887 g CO₂ (United States Environmental Protection Agency (US EPA), 2023).

This study assumes that core refusal in sand marks the Pleistocene surface, and therefore the beginning of the Holocene marsh record. But this interpretation carries some uncertainty. While the basal sand observed at Boundary Bay is likely of Pleistocene origin, it is important to note that other sources of sand, such as tsunami events and fluvial deposition through Fraser River floods, could have influenced the sediment composition of the marsh (Pilarczyk et al., 2021; Dashtgard, 2011). Typically, the true Pleistocene surface lies much deeper in coastal Lower Mainland stratigraphy and is characterised by blue grey glaciomarine clay or glacially derived till (Bednarski, 2015), which is distinguishable from event-related sand deposits caused by tsunamis, storms, rivers, and anthropogenic change (e.g., highway construction). With this information in mind, we note the possibility that the sand layer encountered at the bottom of cores extracted in this study may not represent the true onset of the marsh record.

4.2 Carbon accumulation rates and marsh processes

While the average Boundary Bay marsh CAR (80 ± 45 g C m^-2yr^-1) is lower than CARs measured in nearby Clayoquot Sound, BC (185 ± 50 g C m^-2 yr^-1; Chastain et al., 2022), this average is comparable to both the average global CAR for salt marshes estimated by the Intergovernmental Panel on Climate Change (91 ± 19 g C m^-2 yr^-1; Kennedy et al., 2013), and regional averages of ²¹⁰Pb-dated marshes for the Pacific coast of North America (112 ± 12 g C m^-2 yr^-1, Chastain et al., 2022). We note, however, that Boundary Bay marsh exhibits strong variability in CARs, ranging from 24 ± 12 to 152 ± 51 g C m⁻² yr^-1, with no statistically significant differences between the western and eastern portions of the marsh, largely because of this high spatial variability.

Previous research has suggested that the salt marsh in western Boundary Bay has expanded over the past 90 years (Gailis et al., 2021), but this trend does not appear to extend to the eastern portion of the marsh. Several lines of evidence indicate that the eastern portion of the marsh is eroding, likely in response to a combination of hydrodynamic activities, vegetation, low sediment discharge, and anthropogenic disturbances. For instance, average CARs in BBE (35 ± 11 g C m^-2 yr^-1) are substantially lower than the average for Boundary Bay marsh (80 ± 45 g C m^-2 yr^-1, Table 2). Additionally, the extremely low unsupported ²¹⁰Pb activity in low marsh core BBE1L is consistent with the strong variability within the marsh and supports the hypothesis that the low marsh in eastern Boundary Bay is degrading. Comparisons with aerial photographs further supports this interpretation as the BBE1L coring site was once part of the upper low marsh in 1930 (Figure 6a), whereas in present-day satellite imagery, it lies near the boundary of vegetated marsh (Figure 6b), and currently has a low-density monoculture of Salicornia spp. In fact, Swinbanks and Murray (1981) hypothesized that this decrease in marsh area has been ongoing for at least 4000 years.

This lack of stability or growth in the eastern marsh has several potential contributing factors. Vegetation at BBE1L was dominated by Salicornia spp., whose lower productivity and shallow root systems could reduce C production and trapping (Chastain et al., 2022; Kelleway et al., 2017; Gailis et al., 2021). In combination with hydrodynamic activities such as storms, winds, and erosion, this could lead to progressive soil loss at the marsh edge. Furthermore, sediment supply from the neighbouring Serpentine and Nicomekl rivers is relatively low (Swinbanks and Murray, 1981; Dashtgard, 2011), which has been shown to contribute to low CARs (Peck et al., 2020). Anthropogenic disturbance due to Highway 99 construction adjacent to the marsh in 1942 likely further disrupted sediment delivery and hydrological connectivity, contributing to the observed decrease in marsh area.

No sediment cores were dated from the MB study site due to financial constraints. In 2019, preliminary analyses based on surface elevation tables and marker horizons placed at MB by Ducks Unlimited Canada and Smart Shores Inc suggest some evidence to support the accretion of sediments and increasing mean elevation during a three-to-six-month sampling interval over a period of one year (Christensen and Vadeboncoeur, 2020). The use of marker bed horizons as a means of measuring accumulation rates is relatively inexpensive and simple. However, there are several drawbacks, including the potential disturbance of the layer by bioturbation, hydrological activities (tides) and smearing of markers when coring. Furthermore, this method only has a resolution of the order of ± 1 mm and is limited in temporal resolution for detection of longer-term variations in accumulations (Thomas and Ridd, 2004). Indeed, one of our sediment cores from MB (MB1L) holds essentially no C (peat layer = 2 cm; % C = 0.2%; SCD = 0.002 g C cm^-3), suggesting that very limited carbon is accumulating in the lowest part of the marsh in MB. Therefore, further monitoring over longer time periods is necessary to draw conclusions on current sediment dynamics at the leading edge of eastern Boundary Bay after highway construction.

4.3 Carbon storage

Following the approach of Gailis et al. (2021), we estimated total C storage in two parts of the marsh (BBW and BBE) using three different approaches, to better understand how methodology might influence our estimates of marsh volume. At both locations, we found that C storage estimates produced by the three approaches were statistically similar to each other. All methods assume that average core C stocks and SCDs remain constant throughout the marsh area. Given these assumptions, the traditional method is likely an underestimation of C storage since it only factors in marsh depth at the select locations where sediment cores have been collected rather than providing a comprehensive estimation of marsh depth (at more than 100 locations), and therefore volume, when using the kriging volumetric method. On the other hand, the kriging volumetric method requires the most involved field data collection and may also be an overestimation of C storage if the number of sampled sites is small relative to the study area (Zimmerman et al., 1999). Furthermore, the kriging method uses the depth of refusal (DoR) readings throughout the marsh to derive marsh depth. In Boundary Bay marsh, the DoR is sometimes greater than the depth of the basal peat layer measured in our sediment cores. As such, the simplified volumetric method, which utilizes the depth of the basal peat layer, is likely a more accurate estimation of marsh depth compared to DoR. Based on these considerations, this paper presents only the simplified volumetric C storage estimates.

4.4 Sea level rise and Boundary Bay living dike project

Salt marshes can be threatened by sea level rise if their accretion rates do not exceed projected sea level rise (SLR) or if they do not have the space to migrate inland – a phenomenon known as coastal squeeze (Mcleod et al., 2011; Chmura, 2013; Schuerch et al., 2018). If a marsh is drowned and permanently lost, then the C stored in that marsh will be released into the surrounding coastal waters and atmosphere, ultimately contributing to greenhouse gas emissions (Duarte et al., 2013; Sheehan et al., 2019). Assessing whether a salt marsh will persist into the future is essential to verifying if the site is appropriate for a greenhouse gas offset credits or programs (Chmura, 2013; Macreadie et al., 2019). A regional study conducted across western Canada estimates the relative rate of SLR around Vancouver at 0.6 ± 0.1 mm yr^-1 (Mazzotti et al., 2008), which is slower than the global mean rate at 3.7 mm yr^-1 (IPCC, 2021). This lower rate of SLR might mean that salt marshes in the Vancouver area may face minimal loss in area if sufficient accommodation space for inland migration is provided and sufficient sediment load is available for vertical accretion to continue marsh growth. Using the uncompressed depth and estimated basal ages of the peat layer, the sediment accumulation rates of the marsh average 3.0 ± 1.6 mm y^-1 and range from 1.2 to 5.8 mm yr^-1 (Table 2), which is greater than the regional rate of SLR estimated by Mazzotti et al. (2008). We do note, however, that Boundary Bay marsh is constrained on the landward side by a constructed dike designed to protect the neighbouring farmland, housing, and other infrastructure. This dike limits the possible inland expansion Boundary Bay marsh as an adaptive response to sea level rise.

Nature-based adaptation solutions that maximize inland migration of tidal wetlands may help safeguard wetland persistence with SLR. Inland displacement of coastal flood defences such as dikes, and the designation of nature reserve buffers in upland areas surrounding coastal wetlands are possible solutions (Schuerch et al., 2018). However, given competing land-use constraints within Boundary Bay, BC, the proposed Living Dike project, sponsored by the Cities of Surrey and Delta and the Semiahmoo First Nation, aims instead to raise marsh elevation on the seaward side of the existing dike to offset the projected loss of marsh area due to inundation by SLR. This pilot project began enhancing sediment inputs to the marsh in 2023 and plans to deposit sediment in the low marsh areas over three decades (Readshaw et al., 2018). Although not considered as part of the project engineering plan, this marsh enhancement could also result in build-up of C as a mature plant canopy develops and can trap and store more C, as observed in BBM low marsh areas.

The location of the living dike in eastern Boundary Bay overlaps with the eastern BBE and MB study sites. Low marsh areas in eastern Boundary Bay (BBE, MB) made up 62% of its total marsh area, but only 39% of its total C storage. Although substantial, low marsh C storage in eastern Boundary Bay remains significantly lower compared to high marsh, consistent with findings in western Boundary Bay (Gailis et al., 2021). If this nature-based solution were to consider C sequestration as part of its mitigation goals, then sediment amendment in Boundary Bay would need to be geared toward high rather than low marsh development to maximize the co-benefit of the living dike. To increase the marsh’s ability to vertically accrete, the living dike project will need to ensure that the rates of sediment deposition and marsh development are greater than the local rate of SLR.

4.5 Climate change mitigation and wetland management

Salt marshes provide valuable ecosystem services, including habitat for estuarine wildlife, protection from coastal flooding, and C sequestration (Peck et al., 2020; Chmura, 2013; Duarte et al., 2013; Mcleod et al., 2011). The C sequestration and co-benefits of conservation are important due to their potential inclusion in regional, national, and global climate change adaptation and mitigation strategies (Kauffman et al., 2020). Greenhouse gas emissions from wetland conversion to agricultural land have been reported to be as high as 1,067 - 3,003 Mg CO₂e ha^-1, and land-use change sector is one of the top five sources of current GHG emissions (Kauffman et al., 2020). Since European settlement in the US Pacific Northwest, approximately ~85% of vegetated tidal wetlands have been lost to human land use conversion (Mcleod et al., 2011; Brophy et al., 2019). Due to relatively high C stocks, the blue C ecosystems have been recognized in regional and global climate change mitigation strategies (Kauffman et al., 2020; Donato et al., 2011; Duarte et al., 2013; Mcleod et al., 2011).

In 2021, the government of British Columbia dedicated $27 million to 70 watershed and wetland initiatives through Stronger BC: BC’s Economic Recovery Plan (Government of British Columbia, 2021). This Healthy Watershed initiative focuses on urban, rural, and Indigenous communities highly impacted by COVID-19 pandemic to strengthen natural C sinks through restoration and natural flood management. The Government of British Columbia also introduced a voluntary C market that includes wetland restoration, which allows individuals and corporations to offset their greenhouse gas emissions outside of legally mandated reductions (Sheehan et al., 2019). In addition, British Columbia has legislated greenhouse gas reduction targets of 40% below 2007 by 2030, 60% by 2040, and 80% by 2050 (Climate Change Accountability Act, 2007) and introduced a C tax in 2008 that applies to the purchase and use of fossil fuels (Government of British Columbia, n.d.). Currently, salt marshes are not included in any C offset programs or regulated C markets in British Columbia and Canada. However, knowledge of blue C stocks, accumulation rates, and areas of sequestration would be essential for using blue C ecosystems in climate change mitigation policy.

Our evaluation of blue C behavior in Boundary Bay marsh provides some new insights regarding how blue C could be implemented in climate change mitigation as part of planned nature-based solution initiatives. When nature-based solutions such as the living dike project are implemented for climate adaptation purposes, monitoring of C content would be a relatively simple, natural extension of ongoing monitoring that could help these projects document the potential co-benefit of C storage. However, baseline data are clearly needed to understand sedimentary behaviour in dynamic urban systems such as Boundary Bay, to understand where sediments are likely to accumulate and move. That said, the data collection needed to provide baseline information about C sequestration processes - as in Boundary Bay - have the potential to provide essential information about the long-term dynamics of the system. For example, our work shows that western Boundary Bay has expanded by 20% since 1930 and resulted in an additional 1,549-1,698 Mg of C storage there (Gailis et al., 2021). At the same time, activity in eastern Boundary resulted in a 30% loss in marsh area, which we equate with a loss of 450 Mg of C in this region over the same period. This latter value was calculated using the average of the low marsh C stocks from BBE and MB (25 Mg C ha^-1) multiplied by the area lost (18 ha). Knowledge of porewater salinities help to better understand the potential importance of C loss through methane emissions. BBW and BBM are largely polyhaline and therefore are minimal methane sources. In contrast, BBE and MB are mesohaline marshes that are likely emitting methane, thus reducing the total C uptake from this part of the marsh.

In the case of Boundary Bay, our assessment highlights the spatial and temporal variability in the movement of sediments over the past 90–100 years in response to the interactions between human intervention, circulation dynamics, and vegetative responses that have influenced both C sequestration and sedimentation. These changes highlight the importance of understanding the historical context and future long-term monitoring of any site being considered as a nature-based climate solution, as salt marshes are dynamic systems with the potential to change over decades. Assessing the long-term success of the living dike project will require ongoing monitoring of these sedimentary processes, and incorporating C measurements could be a logical extension of any proposed and ongoing monitoring.

5 Conclusion and next steps

This study provided new information on marsh type, area, volume, C stocks, marsh C storage, and C accumulation rates in the eastern portions of Boundary Bay marsh (BBM, BBE, and MB), which were combined with published research (Gailis et al., 2021) on western Boundary Bay, to develop an improved understanding of marsh processes and blue C storage across all parts of the 222 hectares of the marsh.

Specifically, the total C storage down to the peat layer is estimated to be ~17,360 Mg C, with approximately 85% of the C storage in the western portions of the marsh. The 20% expansion of western Boundary Bay has resulted in a net accumulation of 1,549-1,698 Mg C since 1930 (Gailis et al., 2021). In contrast, the marshes in eastern Boundary Bay contain substantially less C, with an estimated 450 Mg C decrease in C storage due to 30% marsh loss, likely from a combination of erosion and highway construction since 1930. Combined, these changes suggest a net growth of 1100–1250 Mg C since 1930, largely due to marsh expansion in western Boundary Bay.

Carbon accumulation rates average 80 ± 45 g C m^-2 yr^-1 across all of Boundary Bay. While this average is comparable to both regional and global estimates, the CARs show strong spatial variability throughout the marsh, with low values in eastern BB being notable. Additionally, C stocks averaged 71 ± 37 Mg C ha^-for high marsh and 41 ± 36 Mg C ha^-1 for low marsh; these values are comparable to recent global average estimates for salt marsh C stocks that have accumulated to a depth of 30 cm (79.2 ± 38.1 Mg C ha⁻¹, Maxwell et al., 2023), likely because of the young age and shallow nature of the Boundary Bay marsh (average depth of basal peat is 24 cm).

To estimate the future potential of Boundary Bay marsh to sequester C, we have calculated (a) the amount of C that could be accumulated down to a 1-m depth, and (b) the 30-year C stock (to address the short-term mitigative potential of the marsh). These estimates suggest that approximately 7,550 Mg C have accumulated in Boundary Bay over the past 30 years. Additionally, if accumulation continues down to a full 1 meter, it has the potential to hold an additional 38,140 Mg C.

At the base of our sediment cores, we observed a distinctive sand layer, which we interpret as representing the pre-marsh substrate and, therefore, the onset of modern marsh formation. Above this sand layer, we consistently observed silt and clay deposits that likely accumulated during the early stages of marsh development, representing initial sediment trapping and gradual elevation gain prior to substantial peat accumulation. These finer sediments provide a foundation for peat-forming vegetation to establish, and their thickness and composition can inform the timing and dynamics of marsh initiation. The overlying peat layer represents sustained organic matter accumulation characteristic of a mature marsh. Thus, the sand, silt, and clay stratigraphy together provides a clear record of marsh initiation and early development, with the base sand layer marking the earliest establishment of the present marsh ecosystem. This basal sand layer could indicate a change in marsh location and drainage following the establishment of dyke or could have formed during a storm event that knocked out the whole marsh (which would be critical information concerning the long-term stability of a marsh used for a nature-based climate solution). Overall, our study indicates that the Boundary Bay is temporally and spatially variable, highlighting the importance of understanding both the historical and ongoing processes within marsh systems that might be considered for blue C initiatives.

Our study takes initial work at Boundary Bay by Gailis et al. (2021) one step further, but on-going monitoring of blue C sequestration and storage potential would be required for its incorporation in provincial C markets and offset initiatives (Gailis et al., 2021; Emmer et al., 2015; IPCC, 2000). Several additional areas of research could further serve as next steps to form a complete C budget of the Boundary Bay salt marsh ecosystem. These steps include: (1) using gas flux chambers to quantify surface-atmosphere greenhouse gas fluxes, combined with measurements of lateral C export and import via tidal waters, to develop a complete greenhouse gas budget of the marsh (e.g. Poffenbarger et al., 2011; Janousek et al., 2021), (2) installing sulfate, salinity, pH, and temperature loggers for year-round measurements to help improve predictions of greenhouse gas emissions, (3) understanding how local SLR affects inland migration of low marsh zone at Boundary Bay, through improved digital elevation modeling (e,g, using topographic LiDAR).

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

KK: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Visualization, Writing – original draft, Writing – review & editing. HB: Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. MP: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Writing – review & editing. CO: Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. This project was supported by a Natural Sciences and Engineering Research Council (Discovery Grant RGPIN342251) to KEK. Parks Canada provided lab space and materials. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Acknowledgments

We thank M. Gailis for helpful discussions as well as A. Martin, D. Bhullar, S. Kumar, C. Basnayake, M. Duncan, J. Levy, J. Morra, and E. Simpson for their assistance in the field.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2025.1681794/full#supplementary-material

References

Amante C. J. (2018). Estimating coastal digital elevation model uncertainty. J. Coast. Res. 34, 1382–1397. doi: 10.2112/jcoastres-d-17-00211.1

Crossref Full Text | Google Scholar

Anderson R. F., Schiff S. L., and Hesslein R. H. (1987). Determining sediment accumulation and mixing rates using Pb-210, Cs-137, and other tracers: Problems due to postdepositional mobility or coring artifacts. Can. J. Fisheries Aquat. Sci. 44, 231–250. doi: 10.1139/f87-298

Crossref Full Text | Google Scholar

Appleby P. G. (2001). “Chapter 9 Chronostratigraphic techniques in recent sediments,” in Tracking environmental change using lake sediments. Eds. W. M. L. and Smol J. P. (Dordrecht, Netherlands: Kluwer Academic Publishers), 171–203). Basin analysis, coring, and chronological techniques.

Google Scholar

Appleby P. G. and Oldfield F. (1978). The calculation of lead- 210 dates assuming a constant rate of supply of unsupported 210Pb to the sediment. Catena 5, 1–8. doi: 10.1016/S0341-8162(78)80002-2

Crossref Full Text | Google Scholar

Arias-Ortiz A., Masque P., Garcia-Orellana J., Serrano O., Mazarrasa I., Marba N., et al. (2018). Reviews and syntheses: Pb-210-derived sediment and carbon accumulation rates in vegetated coastal ecosystems - setting the record straight. Biogeosciences 15, 6791–6818. doi: 10.5194/bg-15-6791-2018

Crossref Full Text | Google Scholar

Bednarski J. M. (2015). Surficial Geology and Pleistocene stratigraphy from Deep Bay to Nanoose Harbour, Vancouver Island, British Columbia Pacific, Sidney, BC, Canada: Geological Survey of Canada. doi: 10.4095/295609

Crossref Full Text | Google Scholar

Brophy L. S., Greene C. M., Hare V. C., Holycross B., Lanier A., et al (2019). Insights into estuary habitat loss in the western United States using a new method for mapping maximum extent of tidal wetlands. PLOS ONE 14 (8), e0218558. doi: 10.1371/journal.pone.0218558

PubMed Abstract | Crossref Full Text | Google Scholar

Callaway J. C., Borgnis E. L., Turner R. E., and Milan C. S. (2012). Carbon sequestration and sediment accretion in San Francisco Bay tidal wetlands. Estuaries Coasts 35, 1163–1181. doi: 10.1007/s12237-012-9508-9

Crossref Full Text | Google Scholar

Chastain S., Kohfeld K. E., Pellatt M. G., Olid C., and Gailis M. (2022). Quantification of blue carbon in salt marshes of the pacific coast of Canada. Biogeosciences. 19, 5751–5777. doi: 10.5194/bg-2021-157

Crossref Full Text | Google Scholar

Chmura G. L. (2013). What do we need to assess the sustainability of the tidal salt marsh carbon sink? Ocean Coast. Manage. 83, 25–31. doi: 10.1016/j.ocecoaman.2011.09.006

Crossref Full Text | Google Scholar

Chmura G. L., Anisfeld S. C., Cahoon D. R., and Lynch J. C. (2003). Global carbon sequestration in tidal, saline wetland soils. Global Biogeochem. Cycles 17. doi: 10.1029/2002GB001917

Crossref Full Text | Google Scholar

Christensen M. and Vadeboncoeur N. (2020). Chapter 2. Mud Bay Monitoring Report. In: City of Surrey, Prioritizing Infrastructure and Ecosystem Risk from Coastal Processes in Mud Bay (PIER), Final Report MCIP15330, Diamond Head: Vancouver, BC, 133 pp.

Google Scholar

City of Surrey (2021). Mud Bay Foreshore Enhancements. Available online at: https://www.surrey.ca/mud-bay-foreshore-enhancements (Accessed January 10).

Google Scholar

Climate Change Accountability Act. (2007). Statute of British Columbia (SBC) 2007, Chapter 42, s 1. Available online at: https://www.bclaws.gov.bc.ca/civix/document/id/complete/statreg/07042_01 (Accessed November 12, 2025).

Google Scholar

Craft C. B., Seneca E. D., and Broome S. W. (1991). Loss on ignition and Kjeldahl digestion for estimating organic carbon and total nitrogen in estuarine marsh soils: Calibration with dry combustion. Estuaries 14 (2), 175–179. doi: 10.2307/1351691

Crossref Full Text | Google Scholar

Dashtgard S. E. (2011). Linking invertebrate burrow distributions (neoichnology) to physicochemical stresses on a sandy tidal flat: implications for the rock record. Sedimentology 58), 1303–1325. doi: 10.1111/j.1365-3091.2010.01210.x

Crossref Full Text | Google Scholar

Dashtgard S. E. and La Croix A. D. (2015). Chapter 3: Sedimentological trends across the tidal-fluvial transition, Fraser River, Canada: A review and some broader implications. Developments Sedimentology 68), 111–126. Available online at: http://dx.doi.org/10.1016/B978-0-444-63529-7.00005-5.

Google Scholar

Donato D., Kauffman J., Murdiyarso D., Kurnianto S., Stidham M., Kanninen M., et al. (2011). Mangroves among the most carbon-rich forests in the tropics. Nat. Geosci 4 (5), 293–297. doi: 10.1038/ngeo1123

Crossref Full Text | Google Scholar

Drever C. R., Cook-Patton S. C., Akhter F., Badiou P. H., Chmura G. L., Davidson S. J., et al. (2021). Natural climate solutions for Canada. Sci. Adv. 7 (23), eabd6034. doi: 10.1126/sciadv.abd6034

PubMed Abstract | Crossref Full Text | Google Scholar

Duarte C. M., Losada I.J., Hendriks I. E., Mazarrasa I, and Marbà. N (2013). The role of coastal plant communities for climate change mitigation and adaptation. Nat. Climate Change 3, 961–968. doi: 10.1038/nclimate1970

Crossref Full Text | Google Scholar

Durham R. W. and Joshi S. R. (1980). The Pb-210 and Cs-137 profiles in sediment cores from Lakes Matagami and Quevillon, northwest Quebec, Canada. Can. J. Earth Sci. 17, 1746–1750. doi: 10.1139/e80-182

Crossref Full Text | Google Scholar

Eakins J. D. and Morrison R. T. (1978). A new procedure for the determination of lead-210 in lake and marine sediments. Int. J. Appl. Radiat. Isotopes 29, 531–536. doi: 10.1016/0020-708X(78)90161-8

Crossref Full Text | Google Scholar

Emmer I., von Unger M., Needelman B., Crooks S., and Emmett-Mattox S. (2015). Coastal Blue Carbon in Practice: A Manual for Using the VCS Methodology for Tidal Wetland and Seagrass Restoration VM0033. Restore America’s Estuaries and Silvestrum, Arlington, VA, USA, 81 pp.

Google Scholar

Engels S. and Roberts M. C. (2005). The architecture of prograding sandy-gravel beach ridges formed during the last holocene highstand: southwestern british columbia, Canada. J. Sedimentary Res. 75, 1052–1064. doi: 10.2110/jsr.2005.081

Crossref Full Text | Google Scholar

Environment and Climate Change Canada (ECCC) (2021). Nature Smart Climate Solutions Fund (NSCSF). Available online at: https://www.Canada.ca/en/environment-climate-change/services/environmental-funding/programs/nature-smart-climate-solutions-fund.html (Accessed December 12, 2021).

Google Scholar

Gailis M., Kohfeld K. E., Pellatt M. G., and Carlson D. (2021). Quantifying blue carbon for the largest salt marsh in southern British Columbia: implications for regional coastal management. Coast. Eng. J. 63 (3), 275–309. doi: 10.1080/21664250.2021.1894815

Crossref Full Text | Google Scholar

Google Earth Pro. Satellite imagery (2019). Google, Maxar Technologies. Accessed using Google Earth Pro version 7.3.3.7786 (Mountain View, CA, USA: Google).

Google Scholar

Government of British Columbia (2020). Boundary Bay wildlife management area (British Columbia Ministry of Environment). Available online at: https://www2.gov.bc.ca/gov/content/environment/plants-animals-ecosystems/wildlife/wildlife-habitats/conservation-lands/wma/wmas-list/boundary-bay (Accessed May 19, 2021).

Google Scholar

Government of British Columbia (2021). Watershed, wetland projects create jobs, protect environment (British Columbia Ministry of Environment). Available online at: https://news.gov.bc.ca/releases/2021ENV0020-000463 (Accessed May 26, 2021).

Google Scholar

Government of British Columbia (n.d). British Columbia Ministry of Environment. Consultation Backgrounder – Carbon Pricing. Available online at: https://www2.gov.bc.ca/assets/gov/environment/climate-change/action/legislation/carbon-pricing-bg.pdf (Accessed May 28, 2021).

Google Scholar

Granse D., Wanner A., Stock M., Jensen K., and Mueller P. (2024). Plant-sediment interactions decouple inorganic from organic carbon stock development in salt marsh soils. Limnol. Oceanogr. Lett. 9, 469–477. doi: 10.1002/lol2.10382

Crossref Full Text | Google Scholar

Heiri O., Lotter A., and Lemcke G. (2001). Loss on ignition as a method for estimating organic and carbonate content in sediments: reproducibility and comparability of results. J. Paleolimnology 25 (1), 101–110. doi: 10.1023/A:1008119611481

Crossref Full Text | Google Scholar

Howard P. J. A. (1966). The carbon-organic matter factor in various soil types. Oikos 15, 229–236. doi: 10.2307/3565121

Crossref Full Text | Google Scholar

Howard P. J. A. and Howard D. M. (1990). Use of organic carbon and loss-on-ignition to estimate soil organic matter in different soil types and horizons. Biol. Fertil Soils (1990) 9, 306–310. doi: 10.1007/BF00634106

Crossref Full Text | Google Scholar

Howard J., Hoyt S., Isensee K., Pidgeon E., and Telszewski M. (2014). Coastal Blue Carbon: Methods for assessing carbon stocks and emissions factors in mangroves, tidal salt marshes, and seagrass meadows Vol. 184 (Arlington, Virginia, USA: Conservation International, Intergovernmental Oceanographic Commission of UNESCO, International Union for Conservation of Nature).

Google Scholar

IPCC, Intergovernmental Panel on Climate Change (2000). Land use, land-use change, and forestry. Eds. Watson R. T., Noble I. R., Bolin B., Ravindranath N. H., Verardo D. J., and Dokken D. J. (Cambridge, United Kingdom: Cambridge University Press), 375.

Google Scholar

IPCC, Intergovernmental Panel on Climate Change (2021). “Summary for Policymakers,” in Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Eds. Masson-Delmotte V., Zhai P., Pirani A., Conorrs S. L., P an C., Berger S., Caud N., Chen Y., Goldfarb L., Gomis M. I., Huang M., Leitzell K., Lonnoy E., Matthews J. B. R., Maycock T. K., Waterfield T., Yelek i O., Yu R., and Zhou B. (Cambridge University Press, Cambridge, UK).

Google Scholar

Irving A., Connell S., and Russell ,.B. (2011). Restoring Coastal plants to improve global carbon storage: reaping what we sow. PloS One 6, e18311. doi: 10.1371/journal.pone.0018311

PubMed Abstract | Crossref Full Text | Google Scholar

Janousek C., Bailey S., van de Wetering S., Brophy L., Bridgham S., Schultz M., et al. (2021). Early post-restoration recovery of tidal wetland structure and function at the Southern Flow Corridor project, Tillamook Bay, Oregon (Corvallis, OR: Oregon State University, Tillamook Estuaries Partnership, Confederated Tribes of Siletz Indians, Institute for Applied Ecology, and University of Oregon) 219 pp.

Google Scholar

Kauffman J. B., Giovanonni L., Kelly J., Dunstan N., Borde A., Diefenderfer H., et al. (2020). Total ecosystem carbon stocks at the marine-terrestrial interface: Blue carbon of the Pacific Northwest Coast, United States. Glob. Change Biol. 26 (10), 5679–5692. doi: 10.1111/gcb.15248

PubMed Abstract | Crossref Full Text | Google Scholar

Kellerhalls P. and Murray J. W. (1969). Tidal flats at boundary bay, fraser river delta, british columbia. Bull. Can. Petroleum Geology 17, 67–97.

Google Scholar

Kelleway J. J., Saintilan N., Macreadie P. I., Baldock J. A., and Ralph P. J. (2017). Sediment and carbon deposition vary among vegetation assemblages in a coastal salt marsh. Biogeosciences 14, 3763–3779. doi: 10.5194/bg-14-3763-2017

Crossref Full Text | Google Scholar

Kennedy H., Alongi D., Karim A., Chen G., Chmura G., Crooks S., et al. (2013). Chapter 4 Coastal Wetlands. In: Hiraishi T., Krug T., and Tanabe K. (eds.), Supplement to the 2006 IPCC Guidelines for National Greenhouse Gas Inventories: Wetlands. IPCC, Switzerland, 4.1–4.39.

Google Scholar

Krishnaswamy S., Lal D., Martin J. M., and Meybeck M. (1971). Geochronology of lake sediments, Earth Planet. Sc. Lett. 11, 407–414. doi: 10.1016/0012-821X(71)90202-0

Crossref Full Text | Google Scholar

Macreadie P. I., Anton A., Raven J. A., Beaumont N., Connolly R. M., Friess D. A., et al. (2019). The future of Blue Carbon science. Nat. Commun. 10 (1), 3998. doi: 10.1038/s41467-019-11693-w

PubMed Abstract | Crossref Full Text | Google Scholar

Macreadie P. I., Costa M. D. P., Atwood T. B., Friess D. A., Kelleway J. J., Kennedy H., et al. (2021). Blue carbon as a natural climate solution. Nat. Rev. Earth Environ. 2 (12), 826–839. doi: 10.1038/s43017-021-00224-1

Crossref Full Text | Google Scholar

Mathieu G. G., Biscaye P. E., Lupton R. A., and Hammond D. E. (1988). System for measurement of 222Rn at low levels in natural waters. Health Phys. 55, 989–992. doi: 10.1097/00004032-198812000-00015

PubMed Abstract | Crossref Full Text | Google Scholar

Maxwell T. L., Rovai A. S., Adame M. F., Álvarez-Rogel J., Austin W. E. N., Beasy K., et al. (2023). Global dataset of soil organic carbon in tidal marshes. Sci. Data 10 (1), 797. doi: 10.1038/s41597-023-02633-x

PubMed Abstract | Crossref Full Text | Google Scholar

Mazzotti S., Jones C., and Thomson R. E. (2008). Relative and absolute sea level rise in western Canada and northwestern United States from a combined tide gauge-GPS analysis. Journal of Geophysical Research: Oceans 113 (C11). doi: 10.1029/2008JC004835

Crossref Full Text | Google Scholar

Mcleod E., Chmura G. L., Bouillon S., Salm R., BjoÈrk M., Duarte C. M., et al. (2011). A blueprint for blue carbon: toward an improved understanding of the role of vegetated coastal habitats in sequestering CO₂. Front. Ecol. Environ. 9, 552–560. doi: 10.1890/110004

Crossref Full Text | Google Scholar

Morton R. A. and White W. A. (1997). Characteristics of and corrections for core shortening in unconsolidated sediments. J. Coast. Res. 13, 761–769.

Google Scholar

Peck E. K., Wheatcroft R. A., and Brophy L. S. (2020). Controls on sediment accretion and blue carbon burial in tidal saline wetlands: Insights from the Oregon coast, USA. J. Geophysical Research: Biogeosciences 125, e2019JG005464. doi: 10.1029/2019JG005464

Crossref Full Text | Google Scholar

Pendleton L., Donato D. C., Murray B. C., Crooks S., Jenkins W. A., Sifleet S., et al. (2012). Estimating global blue carbon emissions from conversion and degradation of vegetated coastal ecosystems. PloS One 7, e43542. doi: 10.1371/journal.pone.0043542

PubMed Abstract | Crossref Full Text | Google Scholar

Pennington W., Cambray R. S., and Fisher E. M. (1973). Observations on lake sediments using fallout Cs-137 as a tracer. Reprint Nature. 242, 324–326. doi: 10.1038/242324a0

PubMed Abstract | Crossref Full Text | Google Scholar

Pilarczyk J. E., Sawai Y., Namegaya Y., Tamura T., Tanigawa K., Matsumoto D., et al. (2021). A further source of Tokyo earthquakes and Pacific Ocean tsunamis. Nat. Geosci. 14 (10), 796–800. doi: 10.1038/s41561-021-00812-2

Crossref Full Text | Google Scholar

Poffenbarger H. J., Needelman B. A., and Megonigal J.P. (2011). Salinity influence on methane emissions from tidal marshes. Wetlands 31, 831–842. doi: 10.1007/s13157-011-0197-0

Crossref Full Text | Google Scholar

Readshaw J., Wilson J., and Robinson C. (2018). Design Basis for the Living Dike Concept, 35 (Vancouver, BC: West Coast Environmental Law). Available online at: https://wcel.org/sites/default/files/publications/2019-livingdikeconceptbrief-final.pdf (Accessed September 18, 2020).

Google Scholar

Ritchie J. C. and McHenry J. R. (1990). Application of radioactive fallout cesium-137 for measuring soil erosion and sediment accumulation rates and patterns: a review. J. Environ. Qual. 19, 215–233. doi: 10.2134/jeq1990.00472425001900020006x

Crossref Full Text | Google Scholar

RStudio Team (2015). RStudio: Integrated Development Environment for R. RStudio, Inc., Boston, MAhttp://www.rstudio.com/.

Google Scholar

Santisteban J. I., Mediavilla R., López-Pamo E., Dabrio C. J., Zapata M. B.R., García M. J.G., et al. (2004). Loss on ignition: A qualitative or quantitative method for organic matter and carbonate mineral content in sediments? J. Paleolimnology 32 (3), 287–299. doi: 10.1023/B:JOPL.0000042999.30131.5b

Crossref Full Text | Google Scholar

Schuerch M., Spencer T., Temmerman S., Kirwan M. L., Wolff C., Lincke D., et al. (2018). Future response of global coastal wetlands to sea-level rise. Nature 561 (7722), 231–234. doi: 10.1038/s41586-018-0476-5

PubMed Abstract | Crossref Full Text | Google Scholar

Sheehan L., Sherwood E., Moyer R., Radabaugh K., and Simpson S. (2019). Blue carbon: an additional driver for restoring and preserving ecological services of coastal wetlands in tampa bay (Florida, USA). Wetlands. 39 (6), 1317–1328. doi: 10.1007/s13157-019-01137-y

Crossref Full Text | Google Scholar

Shepperd J. E. (1981). Development of a salt marsh on the Fraser River delta at Boundary Bay, British Columbia (Vancouver, BC, Canada: The University of British Columbia). doi: 10.14288/1.0052669

Crossref Full Text | Google Scholar

Swinbanks D. D. and Murray J. W. (1981). Biosedimentological zonation of Boundary Bay tidal flats, Fraser River Delta, British Columbia. Sedimentology 28, 201–237. doi: 10.1111/j.1365-3091.1981.tb01677.x

Crossref Full Text | Google Scholar

Thomas S. and Ridd P. (2004). Review of methods to measure short time scale sediment accumulation. Mar. Geology 207, 95–114. 10.1016/j.margeo.2004.03.011

Crossref Full Text | Google Scholar

United States Environmental Protection Agency (US EPA) (2023). Tailpipe Greenhouse Gas Emissions from a Typical Passenger Vehicle, Office of Transportation and Air Quality, EPA-420-F-23-014. Available online at: https://nepis.epa.gov/Exe/ZyPDF.cgi?Dockey=P1017FP5.pdf (Accessed April 23, 2024).

Google Scholar

Worthington T. A., Spalding M., Landis E., Maxwell T. L., Navarro A., Smart L. S., et al. (2024). The distribution of global tidal marshes from Earth observation data. Global Ecol. Biogeography 33, e13852. doi: 10.1111/geb.13852

Crossref Full Text | Google Scholar

Zimmerman D., Pavlik C., Ruggles A., and Armstrong M. P. (1999). An experimental comparison of ordinary and universal kriging and inverse distance weighting. Math. Geology 31, 375–390. doi: 10.1023/A:1007586507433

Crossref Full Text | Google Scholar