Shoreface erosion counters blue carbon accumulation in transgressive barrier-island systems

Landward migration of coastal ecosystems in response to sea-level rise is altering coastal carbon dynamics. Although such landscapes rapidly accumulate soil carbon, barrier-island migration jeopardizes long-term storage through burial and exposure of organic-rich backbarrier deposits along the lower beach and shoreface. Here, we quantify the carbon flux associated with the seaside erosion of backbarrier lagoon and peat deposits along the Virginia Atlantic Coast. Barrier transgression leads to the release of approximately 26.1 Gg of organic carbon annually. Recent (1994–2017 C.E.) erosion rates exceed annual soil carbon accumulation rates (1984–2020) in adjacent backbarrier ecosystems by approximately 30%. Additionally, shoreface erosion of thick lagoon sediments accounts for >80% of total carbon losses, despite containing lower carbon densities than overlying salt marsh peat. Together, these results emphasize the impermanence of carbon stored in coastal environments and suggest that existing landscape-scale carbon budgets may overstate the magnitude of the coastal carbon sink.


Total Organic Carbon (TOC) and Organic Matter (OM) Values and Uncertainties
We analyzed a subset of sediment subsamples for total organic carbon (TOC) on a Costech Elemental Analyzer, model 4010, coupled to a ThermoFisher DeltaV Isotope Ratio Mass Spectrometer to determine conversion factors for organic matter (OM) to organic carbon (OC) (Supplementary Fig. 2).Prior to analysis, we freeze-dried and powdered samples and removed carbonates by adding 2 drops of 1N HCl and drying overnight at 60°C following methods outlined in refs. 1,2.Average analytic precision (2-σ) for replicate measurements of marsh and lagoon TOC, respectively, were 0.06% and 0.18%.Uncertainties for all OC values were based on 95% confidence intervals and were propagated through to uncertainty estimates for final OC erosion rates.Supplementary Figure 2. Marsh-(green) and lagoon-(brown) specific organic matter (OM) and total organic carbon (TOC) values from sediment core subsamples.Linear regressions (solid lines) were forced through the origin.Green and brown shaded windows demarcate 95% confidence intervals.

Depth-Dependent Shoreface Erosion Rates
Morphologic change along the deep shoreface (that is, above the shoreface toe but below the zone of sediment transport by fair-weather waves) in response to sea-level rise (SLR) or to barrier-island transgression has not been well constrained on observational timescales 3,4 .Because sediment erosion and transport decrease following depth-dependent reductions in wave activity, deeper reaches of the shoreface may respond minimally to, and/or temporally lag, changes occurring in the upper or mid-shoreface 4 .This timescale-dependence is partially due to the increase in depth of closure over longer timescales (multi-decadal to centuries), during which less-frequent but more-powerful storms interact with deeper reaches of the shoreface 3,5 .
However, for the upper and middle shorefaces, morphologic change in response to dynamic sea level and barrier-island movement is expected to occur on decadal to centennial timescales 4 .
Given an estimated regional depth of closure (that is, shoreface toe) of 10-15 m for the Virginia Barrier Islands (VBI) 3,6,7 -which is more than 1.5 m greater than the system-wide average depth of the Holocene lagoon (Supplementary Table 5)-we assume that all marsh and lagoon sediments are fully eroded on the beach and shoreface over the multi-decadal to centennial timescales of interest in this study.Accordingly, we assume a fully erosional transgression, in which all backbarrier deposits are eventually eroded along the shoreface; albeit deeper basal lagoonal sediments are not eroded until exposed at depth, some kilometers offshore.
This assumption is supported by recent bathymetric and seismic mapping along the VBI, which demonstrate that-apart from thick inlet-fill deposits-Holocene sediments thin seaward across the shoreface and are fully eroded within 0.5 km of the modern beach, leaving only a discontinuous thin veneer of sandy lag atop Pleistocene deposits on the shallow shelf 8,9 .
Finally, we note that our approximations of carbon erosion are conservative.Firstly, the thickness of shoreface-eroded Holocene sediments is likely a minimum estimate: for only few islands were we able to collect cores that penetrated to underlying Pleistocene deposits.For all others, we use the (likely incomplete) Holocene lagoon unit thickness that was captured; thus, the actual lagoon thicknesses for some islands may exceed those used in our analyses.Secondly, we note that only the carbon losses from Holocene deposits is considered in our analysis.For example, refs. 7,10report likely erosion and reworking of sandy Pleistocene deposits along the VBI shoreface, re-exposing any OC within these sediments to remineralization and/or export; these are unaccounted for in our budgets, but any erosion of these deeper, relatively carbon-poor deposits would add to shoreface carbon erosion.

Calculation of Shoreline-Change Rates, Shoreline Lengths, and Uncertainties
Shoreline-change rates are based on outputs from the Digital Shoreline Analysis System 11 (DSAS).Island-averaged shoreline change rates are determined from either linear regression rates or end-point-rates of shoreline change for all transects for a given island.For shorelines mapped before 1870, we interpolated shoreline position at year 1870 by using the end-point-rate between survey years directly before and directly after 1870.We used these interpolations as the new starting shoreline position for any time periods beginning with 1870.For time periods considered within this study (that is, 1870-2017, 1870-1942, 1942-1994, 1994-2017), we used the average linear regression rate of island transects only when more than four shoreline positions were available for at least 75% of all island transects; otherwise, we took the average end-point rate between earliest and most recent shoreline-position datapoint for all transects along an island.Each method has drawbacks and advantages 12 .Specifically to our work, the maximum end-point-rate in a time period may be the most accurate way to account for total (non-annualized) lagoon OC erosion, as we assume all transgression results in erosion of lagoon sediments that cannot be recovered during periods of island widening.However, the linear regression rate approach to estimating shoreline change is likely the most accurate method because it accounts for four or more shoreline positions within a given time interval.Therefore, we use linear regression rates when possible, recognizing that the result is a conservative approximation of shoreline change rates and therefore a conservative estimate of annual lagoon OC erosion.Refer to footnotes of Supplementary Table 2 for more information about uncertainty calculations.
Because the lagoon OC erosion term (equation 1 in manuscript; Table 1 in manuscript) depends in part on the shoreline change rate, for time periods during which an island experiences an average negative shoreline-change rates (net growth), the lagoon OC erosion term becomes 0 Gg OC yr -1 .Additionally, we subtracted the equivalent lagoon volume covered during a period characterized by net shoreline regression (island widening) from the volume exposed during the next transgressive period, as a conservative accounting of lagoon previously eroded..Shorelines terminate at the point where they shift ~90° away from the predominant open-ocean orientation.Uncertainties are those introduced during mapping of both shoreline terminuses of a given island, following ref. 12 , with shorelines derived from T-sheets (that is, surveys conducted in 1870 and 1942) introducing 11.7 m uncertainty, those from air photos (that is, surveys conducted in 1994) introducing 5.5 m uncertainty, and those derived from lidar (that is, 2017) introducing 2. *Following ref. 12 , uncertainties calculated as the average of all transect 90% confidence intervals between oldest and youngest shorelines in a given range.
† Uncertainties calculated as the average of all transect EPR uncertainties, which, following ref. 12 are equal to the quadrature addition of mapping uncertainties associated with both shoreline positions (that is, oldest and youngest positions within the time range) divided by the difference in years between shoreline surveys.

Marsh-Exposure Rate Methods and Uncertainties
Marsh-exposure rates, ERmarsh, for specific time periods (Supplementary Table 3) were calculated in three steps.First, using ArcGIS shapefiles from refs. 14,15, we calculated the areal extent of historical marsh exposed on the oceanside of the island at each endpoint time (that is, by 1942, 1994, and 2017).For each time period, the 1870 extent of backbarrier marsh constituted the historical marsh used as the baseline, following the methods of ref. 14 .Because islands in their earliest-mapped positions (surveyed between 1851 and 1888) likely were situated at that time atop some (unmeasurable) area of backbarrier marsh (marsh which would later have been exposed and eroded during subsequent transgression), the second step in estimating ERmarsh involved calculating the percent of the 1942 area of each island that overlaid marsh mapped in the original T-sheets.This percentage was then applied to the island area at the earliest mapped year.For example, if 45% of a given island area in 1942 overlapped with marsh mapped in earlier years, we assumed that 45% of that island was underlain by marsh; we simply applied that same percentage to the original 1800s island area (that is, assumed that 45% of the island area was underlain by marsh at the time of earliest mapping).We then assume that this buried marsh from the 1800s was exposed and eroded in accordance with barrier migration over time.Lastly, to estimate marsh-exposure rates, we interpolated the total area of shoreface-exposed marsh over certain time periods (that is, 1870-2017, 1870-1942, 1942-1994, 1994-2017) to an annual rate.
Because islands were first mapped over a course of 37 years (1851-1888), we used the annual marsh exposure rate to interpolate total areal loss between a single, normalized year (1870) and each subsequent year of interest (that is, 1942, 1994, and 2017).Comparing total losses between 1870 and each subsequent year allowed us to calculate total marsh exposure and annual marsh exposure rates for each time period (reported in the first column under each time interval in Supplementary Table 3).
It is possible portions of marsh mapped in the 1800s may have been eroded in the backbarrier before being buried by the adjacent landward-migrating islands: though still contributing to marsh and blue carbon loss, this component would be unrelated to barrier migration.We therefore include in our uncertainty of marsh-exposure rates a maximum potential amount of backbarrier marsh erosion, based on historic VBI marsh area and marsh loss rates reported in ref. 14 , which computes to an annual VBI marsh erosion rate of 0.12%.Backbarrier marsh erosion is expressed on a per-island basis in Supplementary Table 3 under the headers 'Uncertainty of extent of backbarrier erosion.' Uncertainties associated with mapping are reported in Supplementary Table 3 under the column headers 'Uncertainty of exposed marsh area' and are based on the 4.5% uncertainty associated with mapping the historic marsh area reported in ref. 14  ꝉ The difference between total marsh exposed between 1870-1942 and 1870-1994 (that is, 185,209 m 2 ) was subtracted from the area exposed between 1994-2017, as a conservative accounting of marsh eroded from 1870-1942 and therefore not available for erosion from 1994-2017.

Sediment Core Data and Uncertainties
All cores collected as part of this study were opened, described, photographed, and sampled at a 10-cm resolution (for depths less than approximately 55 cm) in marsh peat and at 20-25 cm resolution through underlying lagoon deposits.We analyzed sample aliquots for organic-matter content (loss-on-ignition; LOI) and continued with the methods outlined in the main text.
We extended our analysis of carbon-loss rates to islands for which we were unable to collect new sediment cores through incorporation of island stratigraphy available in the literature (Supplementary Table 5).Values for 'Average of migrating islands' (Supplementary Table 5) were used to inform the OC erosion rate calculation for Ship Shoal (eq.( 1) in manuscript).
Likewise, Parramore averages were used to inform OC erosion rate values for Hog.Values for which Parramore Island data were not available (for example, OC densities and Holocene marsh thickness) were based on values calculated for 'Average of migrating islands.'Values for 'Average of migrating islands' exclude Parramore and Hog, which largely exhibit narrowing rather than migrating behavior 15 .
With the exception of Met_G which reveals an anomalously thick Holocene unit (interpreted as a relict tidal inlet), cores collected from Assawoman, Metompkin, and Smith Islands as part of this study and others did not penetrate to the Pleistocene.Thus, we use all cores to inform stratigraphy and consider our lagoon thickness estimates for these islands to be conservative.This is in direct contrast to Cedar, where although cores collected as part of this study did not penetrate to the Pleistocene, we can approximate lagoon thickness from the two cores from the literature that did.Sandy facies interbedded within the lagoon unit-likely relict washover fans or tidal inlets-were considered a part of the "lagoon" deposits and thus integrated into the lagoon thickness and average OC density values, increasing the former but decreasing the latter (Fig. 1C).
Shoreline lengths are from ref.
3 m uncertainty.Negative values indicate regression (that is, seaward translation) of the oceanside shoreline (beach growth), while positive values indicate transgression (erosion and/or island migration).
(see Supplementary Table3footnotes for more details).Refer to notes in table above.*A greater area of marsh became exposed seaward of Myrtle from 1870 to 1942 than from 1870 to 1994, leading to a calculation of 0 m 2 marsh loss from 1942 to 1994.