Study area

The Stillaguamish River discharges into the Port Susan Bay Preserve, owned by The Nature Conservancy. The Stillaguamish River drains a 1800 km² watershed largely fed by rain and snowmelt from the North Cascade Range and the Puget lowlands. This mesotidal estuary (2.3 m mean tidal range) faces south-southwest, almost directly facing the direction of predominant winter storm winds, resulting in a storm wave fetch of up to 35–50 km. The nearest tide station provides a local rate of relative sea level rise of 1.89 mm/yr over the period 1972–2019 [14]. The Port Susan Bay Preserve contains a 60-hectare marsh restoration site that was reintroduced to tide and river pulses via dike lowering and two dike breaches in 2012. Prior to the restoration, the site had been managed as diked and drained farmland since the 1950s, and during this period it had subsided up to 1 m in elevation. Since 2011, the wetlands research group at Western Washington University (WWU) has partnered with The Nature Conservancy to conduct pre- and post-restoration monitoring both inside and outside the restoration area (Fig 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Location of (A) Stillaguamish estuary within Washington state, USA, (B) monitoring zones within the estuary, and (C) study sites for this study.

Study sites are denoted with yellow circles (HM = high marsh, LM = low marsh, R = restored, TF = tidal flat). Zone 2 is the 60-hectare restored marsh. ArcGIS World Imagery basemap source: Esri, Maxar, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS user community.

https://doi.org/10.1371/journal.pone.0257244.g001

Beginning in 2011, we established study sites in five monitoring zones (Fig 1). This included the establishment of 21 permanent surface elevation tables (SETs) [15] to measure long-term elevation change throughout the estuary. These zones represent a range of conditions from a healthy, undisturbed marsh (Zone 1), to a regenerating, rapidly accreting emergent marsh in the restoration zone (Zone 2). The study zones also vary hydrologically, with some located directly adjacent to the Stillaguamish River mouth (shown as Hat Slough) and others more distant from the river and its active distributaries (e.g., Zones 4 and 5).

Except for the restoration zone, each of the zones consist of two sites (high and low marsh relative to MHW) and each site has two SETs. The restoration zone consists of four low marsh sites each containing one SET; there is currently no high marsh there due to the 0.5–1.0 m of subsidence prior to restoration. An additional site is located in the lower-elevation tidal flats for a total of 13 sites.

Field sampling

Field sampling permission was provided by The Nature Conservancy as the owner of our entire study area. In the spring and summer of 2016, at each of our 13 sites, we collected two replicate cores to determine soil carbon stocks and vertical accretion rates. At sites containing two established SETs, each replicate core was paired with a SET and located randomly within 10 m of that SET. At sites containing only one SET, replicate cores were both located randomly within 10 m of the same SET. PVC coring tubes (10.2 cm internal diameter) were driven up to one meter into the sediment with a sledgehammer. Cores were kept upright during transport to the laboratory and frozen for later processing.

Aboveground and belowground plant biomass samples were collected at the end of the 2016 growing season (late August) from three replicate plots (0.15 m²) distributed randomly at each site. All aboveground biomass within each plot was clipped to the sediment surface. Belowground biomass samples were collected by pushing PVC coring tubes (10.0-cm internal diameter) to a depth of 50 cm in the center of each replicate biomass plot to capture the entire root zone. Belowground biomass cores were extracted in the field and sectioned into 10-cm depth sections. Biomass samples were stored in plastic bags for transport to the laboratory and temporary storage at 2˚C.

We measured elevation change at the 21 previously established SETs in the spring of each year through 2018, to continue building a time series of elevation change that began in 2011. Field measurement methods followed those described by Lynch et al. [15] for deep RSETs. Elevation change rates were determined using an ordinary least squares linear regression, with time as the independent variable, elevation relative to initial as the dependent variable, and the slope of the regression indicating the rate of surface elevation change. Elevation was previously measured at each SET in 2012 by the US Geological Survey (USGS) using Real Time Kinematic Global Positioning System (RTK-GPS). For this study, where sites included a pair of SETs, we averaged the two SET elevations to report one elevation for each site.

Sediment porewater salinity was measured in late summer of 2016 at the time of plant biomass collection. At each site, one of the pits remaining from the belowground biomass cores was left open to allow porewater to collect. Porewater salinity was measured in the field with a portable refractometer. A few sites were completely dry during sampling, preventing salinity measurements at those sites.

Laboratory analyses

Plant aboveground biomass samples were washed with a 0.5-mm sieve and separated by species. Belowground biomass samples were washed and separated into live and dead root fractions, with the dead fraction then discarded to isolate end-of-season live biomass. Dead roots present in various states of decay and are thus considered a component of the soil fraction. All samples were oven dried at 60 ˚C for at least 96 hours and weighed. We used a 1:0.38 ratio to convert biomass dry weight to biomass carbon [16].

In the laboratory, frozen sediment cores were extracted, sliced into 2-cm sections with a power saw, and dried at 60 ˚C for at least 96 hours to determine bulk density. A subsample of each section was pulverized using a Wiley Mill with 0.425-mm mesh screen. The organic matter (OM) content of each section was determined by loss on ignition (LOI). Subsamples were burned at 500 ˚C for approximately 8 hours and weighed before and after burning [17].

We measured percent organic carbon (C_org) content for a subset of samples (n = 124) with a FlashEA 1112 CN analyzer (Thermo Fisher Scientific, Waltham, MA). We first measured total (organic and inorganic) carbon content, then measured inorganic carbon by analyzing the ashed subsamples that remained after LOI, and finally calculated the C_org content by subtraction. We developed an OM-C_org conversion from these 124 samples that was then applied to all sediment samples to produce the C_org contents reported here (C_org = 0.0035*OM² + 0.4135*OM– 0.4496; R² = 0.97; P < 0.001) (S1 Fig). We used C_org for all further calculations of carbon values (carbon density, carbon stocks, carbon accumulation rates). Carbon density was calculated for each core section as the product of C_org content and bulk density. Each sediment characteristic parameter was averaged across the top 30 cm of each core for reporting purposes. Sediment carbon stocks were also calculated to a depth of 30 cm. However, a thick layer of organic material (presumably from a depositional event or burial of pre-restoration biomass) was observed in the restoration site cores between 10 cm and 20 cm in depth, which biased those results high when averaged across the top 30 cm. We therefore used only the top 10 cm of the restoration site profiles and extrapolated to 30 cm to best represent post-restoration conditions without the influence of this depositional layer.

We analyzed sediment grain size composition by separating fine from coarse sediments to evaluate the relationship of grain size with sediment carbon stocks. A subset of six samples distributed in even depth intervals was selected for analysis from each core. Each 30-mL subsample was redried, weighed, soaked, and rinsed through a 63-μm sieve to remove the fine sediment (silt and clay) portion. The remaining coarse-grained (sand) portion was dried and weighed.

Long-term sediment accretion rates were determined from the downcore distribution of excess ²¹⁰Pb activity from one core per site. Our analysis was limited to one core per site due to time constraints (the analysis can require up to one month per core). We used the Canberra Germanium Detector (model GL2820R, Mirion Technologies (Canberra) Inc., Meriden, CT), with gamma emissions at 46 keV and 351 keV recorded by Genie 2000 software (Canberra 2002). Excess (unsupported) ²¹⁰Pb was calculated as the difference between total ²¹⁰Pb (at 46 keV) and supported ²¹⁰Pb (at 351 keV) to distinguish between excess ²¹⁰Pb deposited at the sediment surface and supported ²¹⁰Pb that has decayed from radium in the sediment. Within each core, a subset of 2-cm sections was analyzed from the surface to the depth at which excess ²¹⁰Pb declined to zero. We used the constant initial concentration (CIC) model [18] because in many cases we could not recover the full profile of excess ²¹⁰Pb due to core length or time constraints. With CIC, a linear regression of the natural log of excess ²¹⁰Pb activity versus depth was used to determine the sediment accretion rate, which is equal to –λ/s, where λ is the half-life of ²¹⁰Pb (22.2 yr^-1) and s is the slope of the regression. Because ²¹⁰Pb profiles from the four restored sites showed evidence of disturbance, preventing us from assigning them a long-term accretion rate that could be compared with rates from the undisturbed sites, we used the mean ratio of ²¹⁰Pb-based rate to SET-based rate from the eight natural marsh sites to estimate the ²¹⁰Pb-based accretion rate for the restoration sites.

Carbon accumulation rates were calculated in two ways: 1) as the product of sediment carbon density and the ²¹⁰Pb -based accretion rate, and 2) as the product of carbon density and the SET-based elevation change rate. These two rates are reported as the ²¹⁰Pb-based carbon accumulation rate and the SET-based carbon accumulation rate, respectively. Although SET-based carbon accumulation rates are not as commonly reported in blue carbon studies, we hypothesized that, in concept, these rates should be comparable to ²¹⁰Pb -based accumulation rates. We include them to 1) provide a validation for the ²¹⁰Pb -based rates, 2) test the feasibility of using SET and ²¹⁰Pb accretion rates interchangeably, and 3) provide a potential alternative to the ²¹⁰Pb -based rates in the event that any cores were unable to be dated with ²¹⁰Pb due to sediment mixing or other disturbance.

Statistical analyses were conducted in the R programming environment using a critical value (α) of 0.05 [19]. We used the Shapiro-Wilk test to check for normality of biomass and core data. We used simple linear regression to determine the relationship between organic content and C_org content. We used Spearman’s rank correlation to evaluate relationships between environmental variables and plant biomass, sediment carbon stocks, and carbon accumulation rates, and Kruskal-Wallis H test to compare high and low marsh biomass.

Plant biomass and biomass C stocks

Aboveground biomass ranged from 0 g dw m^-2 (at Site TF, the tideflat site) to 978.64 g dw m^-2 (HM5) (Fig 2). Aboveground biomass at low marsh sites averaged 279.09 g dw m^-2, compared to 879.40 g dw m^-2 at high marsh sites, with 353.37 g dw m^-2 across the four restoration area sites. Belowground biomass across all sites was an average of 2.4 times greater than aboveground biomass (median of 1.6), ranging from 0 g dw m^-2 (TF) to 2,243.98 g dw m^-2 (HM5). Belowground biomass at low marsh sites averaged 397.83 g dw m^-2, compared to 1,686.54 g dw m^-2 at high marsh sites, with 329.11 g dw m^-2 at restoration area sites.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Aboveground and belowground biomass by site.

Sites are arranged in order of increasing elevation. The elevations of the restored sites (R1 –R4) were within the range of other low marsh sites. Error bars represent standard error.

https://doi.org/10.1371/journal.pone.0257244.g002

Sediment characteristics

The mean sediment bulk density observed at both natural and restored marsh sites was identical, at 0.92 g cm^-3, lower than the tidal flat bulk density of 1.48 g cm^-3 (Table 1). In general however, organic matter content, organic carbon content, and carbon density values were highest at natural marsh sites, intermediate at restored sites, and lowest at the tidal flat site. Organic matter content averaged 6.45% at natural marsh sites, 4.93% at restored sites, and 2.70% at the tidal flat site. Carbon content averaged 2.38% at natural marsh sites, 1.74% at restored sites, and 0.70% at the tidal flat site. Carbon density averaged 0.020 g C cm^-3, 0.015 g C cm^-3, and 0.010 g C cm^-3, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Site elevation, marsh type, porewater salinity, and sediment characteristics.

Sediment characteristics (bulk density, organic content, carbon content, carbon density, and grain size) are reported as the average value across the top 30 cm. Elevation is given relative to mean lower low water (MLLW).

https://doi.org/10.1371/journal.pone.0257244.t001

The percentage of fine-grained sediments, as averaged over the top 30 cm, was lowest (23.46%) at the tidal flat site (TF), and highest (98.60%) at high marsh site HM4. Although the percentage of fines was more variable across the natural marsh sites, the mean from both natural marsh sites and restored sites was approximately equal, with 76.02% and 76.66% fines, respectively (Table 1). As expected, sediment grain size from individual samples was significantly correlated with sample organic content (Spearman’s rank correlation, ρ = 0.78, P < 0.001, n = 151). However, a visual assessment of the relationship suggests that the percentage of fine-grained sediments best correlated with the minimum organic content rather than the full range of observed values at any given site (S2 Fig).

Elevation change rates

Rates of surface elevation change measured with SETs ranged from -0.03 cm yr^-1 (TF) to 3.79 cm yr^-1 (R2) (Table 2). The mean rate across the eight natural marsh sites was 1.03 cm yr^-1 (Fig 3), while the mean rate across the four restored sites was nearly three times higher at 2.74 cm yr^-1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Mean rate of elevation change as measured by SETs and averaged by zone, for the entire period of record (2011–2018) with standard error bars.

Zone 2 is the restored area, Zone TF is the tidal flat zone, and the remaining Zones 1, 3, 4, and 5 are natural marsh zones (see Fig 1 for a map showing the location of all zones).

https://doi.org/10.1371/journal.pone.0257244.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Rates of elevation change, accretion, and carbon accumulation, and carbon stocks across the three carbon pools.

SET elevation change rates represent a 4–7 year time period of record, while ²¹⁰Pb accretion rates represent an approximately 100-year record. Both C accumulation rates incorporate the 30-cm average soil carbon density.

https://doi.org/10.1371/journal.pone.0257244.t002

Long-term accretion rates

Long-term accretion rates determined from ²¹⁰Pb dating ranged from 0 cm yr^-1 (TF) to 2.17 cm yr^-1 (R2) (Table 2). The mean accretion rate at natural marsh sites was 0.60 cm yr^-1, while the mean rate at restored sites was 1.57 cm yr^-1. The restoration site ²¹⁰Pb profiles appeared disturbed, with discontinuous layers and a mixed surface layer corresponding to the most recent post-restoration period, therefore the mean ²¹⁰Pb:SET ratio of 0.57 was applied to SET-based rates at these four sites to estimate the ²¹⁰Pb accretion rate. ²¹⁰Pb accretion rates varied more among zones than within zones, although a Kruskal-Wallis test determined this pattern to be not significant (H = 10.14, 5 df, P = 0.07).

Most cores demonstrated anomalous ¹³⁷Cs profiles with low activities and no obvious peak, indicating either post-depositional vertical mobility [20] or peak deterioration due to radionuclide decay [21], therefore we do not report any ¹³⁷Cs-derived accretion rates.

Carbon stocks and accumulation rates

Sediment carbon stocks in the top 30 cm ranged from 2.95 kg C m^-2 (TF) to 7.95 kg C m^-2 (HM5) (Table 2). As with other sediment characteristics, restored site carbon stocks were intermediate between the tidal flat site and the natural marsh sites, averaging 4.43 kg C m^-2 in restored sites, 2.95 kg C m^-2 in the tidal flat, and 5.95 kg C m^-2 in natural marshes. Sediment carbon stocks were approximately 30 times greater on average than biomass carbon stocks (Fig 4, Table 2). Aboveground biomass carbon stocks ranged from 0 kg C m^-2 (TF) to 0.37 kg C m^-2 (HM5 and R4), averaging 0.22 kg C m^-2 at the natural marsh sites, and 0.13 kg C m^-2 at the restored sites. Belowground biomass carbon stocks ranged from 0 kg C m^-2 (TF) to 0.85 kg C m^-2 (HM5), averaging 0.40 kg C m^-2 at natural marsh sites and 0.13 kg C m^-2 at restored sites.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Carbon stocks (±SE) in aboveground biomass, belowground biomass, and the top 30 cm of sediment at each site.

Sites are arranged in order of increasing elevation.

https://doi.org/10.1371/journal.pone.0257244.g004

Carbon accumulation rates calculated with ²¹⁰Pb-based accretion rates ranged from 0 g C m^-2 yr^-1 (TF) to 350.23 g C m^-2 yr^-1 (R2), averaging 123.00 g C m^-2 yr^-1 at natural marsh sites and 230.49 g C m^-2 yr^-1 at restored sites. Carbon accumulation rates calculated with SET-based elevation change rates ranged from -2.95 g C m^-2 yr^-1 (TF) to 611.69 g C m^-2 yr^-1 (R2), with an average of 207.82 g C m^-2 yr^-1 at natural marsh sites and 402.93 g C m^-2 yr^-1 at restored sites (Fig 5, Table 2). Although the ²¹⁰Pb and SET-based rates represent different time periods of record (100 years and 4–7 years, respectively), the same 30-cm depth was used for mean carbon density in both cases.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Carbon accumulation rates (±SE) from each site, with both ²¹⁰Pb-based and SET-based rates shown.

Sites are arranged in order of increasing elevation.

https://doi.org/10.1371/journal.pone.0257244.g005

Environmental characteristics

Porewater salinity, as measured at the end of the growing season in mid-August, ranged from 6 ppt (at HM1 adjacent to Hat Slough) to 20 ppt (at several sites), and averaged 15.5 ppt within both natural marsh sites and restored sites (Table 1).

Site elevations ranged from 2.15 m MLLW (TF) to 3.34 m MLLW (HM5). Undisturbed low marsh sites ranged from 2.48 to 2.71 m MLLW, while high marsh sites were approximately 0.5 m higher, ranging from 3.04 to 3.34 m MLLW. Elevations within the restoration area fell within the range of undisturbed low marsh elevations, ranging from 2.40 to 2.85 m MLLW (Table 1).

Drivers of plant biomass, carbon stocks, and carbon accumulation rates

Plant biomass was more closely related to elevation than to salinity. Elevation positively correlated with aboveground biomass (Spearman’s rank correlation, ρ = 0.86, P = 0.0002, n = 13) while salinity did not correlate with aboveground biomass (Spearman’s rank correlation, ρ = -0.27, P = 0.45, n = 10). Among the natural marsh sites, Kruskal-Wallis tests determined that both aboveground and belowground biomass were significantly greater at high marsh sites compared to low marsh sites (Aboveground: H = 5.33, 1 df, P = 0.02; Belowground: H = 5.33, 1 df, P = 0.02).

Sediment carbon stocks were significantly and positively correlated with elevation, aboveground biomass, and belowground biomass, but the strongest relationship was with elevation (Spearman’s rank correlation, ρ = 0.86, P = 0.0003, n = 13). Carbon accumulation rates, however, only correlated with accretion rates such that high accretion rates led to high carbon accumulation rates (Spearman’s rank correlation, ρ = 0.97, P < 0.001, n = 13). These two rates are not independent however, with carbon accumulation calculated as the product of accretion rate and carbon density. No significant relationship was detected between carbon accumulation rates (both SET-based and ²¹⁰Pb-based) and elevation, biomass, or even sediment carbon stocks (which are based on carbon density).

Plant biomass

Within natural marshes, aboveground biomass at high marsh sites was approximately three times greater than at low marsh sites. Belowground biomass showed even more of a contrast, with over four times greater biomass at high marsh sites. The restored sites varied considerably in terms of both aboveground and belowground biomass. Restored site R4 had a higher aboveground biomass than other low marsh sites. This was the southernmost site in the restored marsh where vegetation was dominated by cattails (Typha latifolia and T. angustifolia) and maritime bulrush (Bolboschoenus maritimus), both relatively productive species. Cattails were not found at any other site. Their prevalence in the restoration area is likely a result of their occurrence across the site’s freshwater wetlands prior to restoration, and their ability to persist after restoration despite suboptimal conditions [22]. Restored sites R1 and R2 contained very little aboveground biomass. These two sites were located near the northern dike breach where a developing sand splay reduced the area’s capacity for permanent vegetation establishment.

Sediment characteristics

Sediment profiles at all four restored sites showed a spike in organic material between 15 and 20 cm in depth. The cause of this spike is unclear but it is hypothesized to be from either pre-restoration standing plant material that was quickly buried after the dike breach (or trampled as part of the levee removal process), or a layer of material originating outside the site that was deposited during a storm event. The latter is perhaps more likely since this layer looked more like small pieces of wood instead of dead herbaceous material. In any case, the depth of this layer reasonably corresponds to the depth of sediment we expect to have accumulated after the dike breach, based on SET elevation change rates from the last four years.

All other undisturbed sites (beyond the restoration area) generally exhibited typical sediment profile patterns with highest percent organic matter and carbon at the surface, and a decrease with depth due to organic material decomposition and diminishing root growth. Low marsh and tidal flat sites tended to have lower levels of carbon at the surface, but they still showed some decrease with depth. Bulk density was generally lowest at the surface as expected, increasing with depth due to the compaction and decomposition of less dense organic matter.

The restored marsh sites, although low in elevation, had a mean sediment bulk density comparable to that of the high marsh sites. Despite the restoration and natural marsh sites having a similar mean bulk density and percentage of fine-grained sediments, organic content, carbon content, and carbon density were all lower in the restoration area. In fact, these values were low in the restored sites even when compared to just the natural low marsh sites, which occur within the same elevation range (e.g., restored site carbon density averaged 0.015 g C cm^-3 compared to 0.017 g C cm^-3 in natural low marshes), Because the restoration area has experienced very rapid sediment deposition since the dike breaches in 2012, plant material has had little opportunity to develop and influence the sediment profile. However, even the sites with little or no plant material (TF, R1, and R2) still contained a small amount of sediment carbon, indicating that this carbon is likely allochthonous, presumably deposited with the incoming sediment independent of the in situ vegetation. Similarly, Wollenberg et al. [23] measured a 2.3% carbon content in newly deposited sediment in the Bay of Fundy, concluding that short-term soil carbon density was driven by the inherent carbon content of incoming sediment.

Sediment carbon density in the Stillaguamish River estuary is low relative to many other estuaries; indeed the average values here are lower than the United States mean of 0.027 g cm^-3 [24]. Sediment organic matter content, carbon content, and carbon density across the Stillaguamish River estuary are also slightly lower than those reported for the neighboring Snohomish River estuary near Everett, WA [12]. For example, carbon density at six Snohomish tidal marsh sites ranged from 0.018 g C cm^-3 to 0.032 g C cm^-3, whereas Stillaguamish marsh carbon density ranged from 0.012 g C cm^-3 to 0.026 g C cm^-3. Oregon tidal marsh carbon density is higher than both these Washington estuaries, with an Oregon mean of 0.034 g C cm^-3 [25] compared to this study’s natural marsh mean of 0.20 g C cm^-3. Although these two neighboring estuaries are only separated by approximately 20 miles, the most likely explanation for the difference in sediment characteristics is their differing hydrodynamics. The Stillaguamish estuary marshes are located at the northern end of Port Susan Bay with full exposure to waves from winter storms. In contrast, the Snohomish estuary marshes are situated further inland, where they are more protected from storms, allowing finer sediments to settle. Having blue carbon data from both estuaries provides a window into the variability we might expect to see across the Salish Sea, in terms of both storm exposure and access to riverine sediments.

Accretion and elevation change rates

All sites, with the exception of the tidal flat site, are accumulating sediments and gaining elevation at a rate that exceeds current and even some predicted future rates of sea level rise, indicating that the estuary is not currently sediment-limited.

SETs in Zones 1 and 3 have seen a fairly steady increase in elevation over their 7-year period of record. These zones are nearest to the Stillaguamish River mouth and riverine sediment supply. Previously at Zones 4 and 5 (farthest from the river), sediment starvation and wave-induced erosion were believed to be deteriorating the low marsh and causing the high marsh to retreat. However, these areas have seen very dynamic conditions in the most recent years, with pools and channels encroaching in one year, and an influx of large woody debris and sediment deposits in the next. Although these dynamic conditions lead to more intermittent deposition rather than steady and predictable accretion, they appear to be slowing erosion in these areas.

The Zone 2 restored site SETs have shown 2.74 cm yr^-1 of elevation change on average, over their four-year period of record. This rate is nearly three times that of the other undisturbed marsh sites. While it was expected and desired that the return of tidal and riverine pulsing to this formerly hydrologically isolated and subsided farmland would result in enhanced sedimentation, the rates observed here are remarkably high. These rates are higher than the 0.5–0.8 cm yr^-1 calculated by Nowacki and Grossman [26] based on their measured suspended sediment fluxes and our mean bulk density, suggesting that perhaps an assumption of uniform sediment distribution is not accurate. Despite our initial high rates of elevation change, it is expected that as infill continues and the elevation inside the restoration zone approaches that of the adjacent undisturbed marshes, rates of elevation change will gradually decline. These initial high rates demonstrate that the restored marsh has the capacity to restore its elevation without additional human intervention, and this is likely to apply to other Puget Sound estuaries with undammed mountainous watersheds.

We originally predicted that SET-based elevation change rates and ²¹⁰Pb-based accretion rates would be comparable because they measure similar processes of sediment deposition and shallow subsidence. However, we found the ²¹⁰Pb accretion rates to be consistently lower than SET rates, with an average ²¹⁰Pb:SET ratio of 0.57 across the eight natural marsh sites. One possible explanation for this difference is that the ²¹⁰Pb method does not measure the most recent depositional layer with as much precision as the SET method, providing only one data point over the top 2 cm. This may result in the most recently deposited sediment not being fully accounted for, which would reduce the reported rate. Additionally, this surface layer may be compacted somewhat during the coring process. Another possible reason for the discrepancy is that SETs have a shorter length of record. Very few studies have used both methods at the same sites to allow for a comparison, but Breithaupt et al. [27] also found that SETs recorded a higher rate of accretion than that measured by ²¹⁰Pb. Our results indicate that the two methods cannot be used interchangeably in the Stillaguamish estuary, although one rate may potentially be predicted from the other with this site-specific ratio.

Carbon stocks and accumulation rates

Global average rates of salt marsh carbon accumulation are generally higher than those reported here for Stillaguamish natural marshes. For example, Mcleod et al. [2] reported a global average rate of 218 g C m^-2 yr^-1. A more recent review by [28] reported a global average of 244.7 g C m^-2 yr^-1, but a lower NE Pacific regional average (based on eight California sites) of 173.6 g C m^-2 yr^-1. However, some of the studies included in those global reviews used shorter-term accretion methods such as feldspar, which overestimate long-term accretion. Other studies along the west coast of North America report natural marsh rates more comparable to ours (all using ²¹⁰Pb dating), including an estimated 79 g C m^-2 yr^-1 from San Francisco Bay [29], 77 g C m^-2 yr^-1 from Oregon [25], and 115 g C m^-2 yr^-1 from British Columbia, Canada [30]. Within the Puget Sound, Crooks et al. [12] reported 110 g C m^-2 yr^-1 from a natural marsh in the Snohomish estuary and Drexler et al. [13] reported a 50-year rate of 118 g C m^-2 yr^-1 from the Nisqually River estuary, both similar to our Stillaguamish average of 123 g C m^-2 yr^-1.

Few studies report carbon sequestration rates following marsh restoration. In Oregon, Brophy et al. reported 80 g C m^-2 yr^-1 in Columbia River estuary restored sites [10] and 150 g C m^-2 yr^-1 in Tillamook Bay restored sites [11], both lower than Stillaguamish rates. Drexler et al. [13] reported 81 to 210 g C m^-2 yr^-1 from the Nisqually River estuary. Snohomish estuary rates, however, were most similar to Stillaguamish rates, with up to 352 g C m^-2 yr^-1 based on 1.61 cm yr^-1 of accretion [12]. Of course, it is not fair to compare one restored site to another without accounting for differences in vertical accommodation space, which is affected by site subsidence prior to restoration and time since restoration, but such an analysis is beyond the scope of this study.

Carbon stocks and accumulation rates are driven by different processes in the Stillaguamish River estuary. In other words, a high carbon stock does not necessarily imply a high carbon accumulation rate. Carbon stocks were highest at high marsh sites, whereas carbon accumulation rates at those same sites were moderate. Stocks were most strongly related to elevation (as a proxy for tidal inundation frequency), whereas accumulation rates were only predicted by sediment accretion rates, similar to the findings of Wollenberg et al. [23] and Peck et al. [25]. Accretion rates themselves are quite variable throughout the Stillaguamish estuary, ranging from 0.24 to 1.24 cm yr^-1 just within natural marsh sites. These rates are clearly driven not just by estuary-scale conditions but also by smaller-scale sediment dynamics that warrant further study.

When comparing restored sites to just the natural low marsh sites, we might expect those carbon stocks and accumulation rates to be more similar because those sites exist within the same elevation range. However the mean sediment carbon stock in restored sites (4.43 g C m^-2) was still lower than in the natural low marsh sites (5.18 g C m^-2), suggesting that sediment carbon density in the restored marsh is still in the process of recovering after restoration. Since the restored marsh has had such high accretion rates from mineral sediment inputs (1.57 cm yr^-1, versus 0.58 cm yr^-1 in natural low marshes), we expect it may be several more years before the restoration soil profile resembles the natural sites.

A few studies have reported carbon accumulation rates calculated with SET elevation change rates [e.g., 31–33] but most use radiometric (²¹⁰Pb or ¹³⁷Cs) accretion rates. We have reported both SET and ²¹⁰Pb rates here to allow for comparisons with both types of studies. SET-based carbon accumulation rates here were all higher, and on average nearly twice the ²¹⁰Pb-based carbon accumulation rates, so they cannot be considered interchangeable. Considering that most studies use ²¹⁰Pb-based carbon accumulation rates, the rest of our discussion emphasizes those rates.

Implications for carbon finance and restoration

Since the Stillaguamish River estuary restored marsh was 0.5 to 1.0 m below the elevation of adjacent natural marshes at the time of dike removal, the restored marsh is expected to continue accumulating sediment until it reaches its former elevation. Assuming a constant sediment carbon density (likely conservative since it is currently lower than adjacent marshes), the entire 60-hectare restored marsh is expected to accumulate approximately 4,500 to 9,000 Tonnes of carbon, which is equivalent to removing approximately 3,500 to 7,000 cars from the roads for one year. That amount of carbon is worth approximately $65,000 to $165,000 using a national average carbon offset price of $4 to $5 per Tonne CO₂. Once its elevation reaches an equilibrium with the surrounding marsh, the restored marsh is expected to continue accumulating carbon at a rate similar to that of the surrounding marsh, estimated at 1.23 Tonnes C ha^-1 yr^-1, or approximately 75 Tonnes C yr^-1 within the entire restoration site. This carbon accumulation rate is valued at approximately $20 ha^-1 yr^-1, or $1,200 yr^-1 within the entire restored area.

The national carbon offset value is currently quite low compared to the estimated social cost of carbon of $41 per Tonne [4]. There is some indication that the carbon offset value awarded to blue carbon projects may increase, in consideration of the many other co-benefits that marsh restoration provides, including but not limited to coastal storm protection, biodiversity, tourism, and fisheries enhancement. Until then, carbon credits alone are unlikely to fund an entire restoration project. However, although the monetary value of carbon sequestration is currently low, there are other potential benefits to recognizing and quantifying the climate benefits of marsh restoration: (1) the climate benefits of restoration may garner interest from additional funders seeking to support a project with a climate change mitigation component; (2) annual carbon offset payments may help fund monitoring efforts, which are typically difficult to fund; and (3) blue carbon adds to a long list of other marsh restoration benefits, increasing justification for restoration.

Tidal wetland restoration projects tend to be relatively small, losing cost-effectiveness once the carbon quantification, feasibility, design, and monitoring are accounted for [e.g., 34]. Fortunately, it is possible to group smaller carbon finance projects together within a larger project area to increase cost-effectiveness. Larger project areas of at least 500 ha could be at the scale of an estuary or even a region, and the individual restorations need not be implemented at the same time. This possibility should be an important consideration for tidal wetland restoration projects.

Recommended next steps

This study expands the amount of PNW blue carbon data available to scientists, managers, and policymakers. To our knowledge no other blue carbon study in the PNW has included carbon accumulation rates from both ²¹⁰Pb and SET measurements, in addition to C stocks within all three marsh carbon pools (i.e., sediment, aboveground biomass, and belowground biomass). The extensive dataset generated by this study can contribute to a feasibility analysis that determines if the site would be appropriate for the carbon market, using the Port Susan Bay Preserve restored marsh as an example for other potential projects within the region. This study will also facilitate comparisons of carbon stocks and accumulation rates across estuaries, which until only recently were not possible due to a lack of data within the region.

Some further research is recommended that was outside the scope of this study, to allow for a full feasibility analysis of the estuary’s carbon finance potential. Greenhouse gas (CO₂, CH₄, and N₂0) emissions must be quantified for a complete blue carbon assessment. Even with substantial carbon sequestration, the estuary’s brackish marshes could simultaneously be emitting greenhouse gases, and it is yet unknown whether the estuary is a net greenhouse gas sink or source. It is possible to use published emission data from other estuaries but published methane (CH₄) emissions can be quite variable, particularly from brackish marshes such as those in the Stillaguamish estuary, therefore site-specific data are preferred.

Carbon offset credits are awarded based on net autochthonous carbon accumulation (carbon produced in situ). Allochthonous carbon (carbon produced outside the system) must be subtracted from the carbon balance for crediting purposes. Further study using stable isotopes would allow for a differentiation between autochthonous and allochthonous carbon to provide this information. We expect that the majority of carbon sequestered in the Stillaguamish estuary restored marsh is allochthonous, given that sediment carbon content is relatively low [35] and a similar restored marsh in Puget Sound reported high allochthonous carbon contribution based on stable isotopes [36].

While much of the sequestered carbon in the restored marsh may be imported from elsewhere, the marsh may also be exporting carbon from its aboveground plant production. The quantity and fate of exported carbon can be difficult to track [37] and it is currently unaccounted for in this study. However, exported carbon can be a substantial contributor to the carbon budget of the larger coastal ocean [38].

Both the Stillaguamish and Snohomish estuaries have demonstrated relatively high rates of accretion (and consequently carbon accumulation) immediately following restoration, higher than rates reported from Oregon estuaries [10, 11] and southern Puget Sound [13] but comparable to rates from the Columbia River estuary [39]. It would be helpful to expand this assessment to restored marshes in other PNW estuaries to determine whether this high capacity for accretion is estuary-specific or region-specific.