Seagrass blue carbon spatial patterns at the meadow-scale

Most information on seagrass carbon burial derives from point measurements, which are sometimes scaled by meadow area to estimate carbon stocks; however, sediment organic carbon (Corg) concentrations may vary with distance from the meadow edge, resulting in spatial gradients that affect the accuracy of stock estimates. We mapped sediment Corg concentrations throughout a large (6 km2) restored seagrass meadow to determine whether Corg distribution patterns exist at different spatial scales. The meadow originated from ≤1-acre plots seeded between 2001 and 2004, so we expected Corg to vary spatially according to the known meadow age at sample sites and with proximity to the meadow edge. Applying spatial autoregressive models allowed us to control for spatial autocorrelation and quantify the relative effects of edge proximity and age on Corg concentrations. We found that edge proximity, not age, significantly predicted the meadow-scale Corg distribution. We also evaluated relationships between Corg and a variety of specific explanatory variables, including site relative exposure, shoot density, sediment grain size, and bathymetry. Factors known to affect carbon burial at the plot-scale, such as meadow age and shoot density, were not significant controls on the meadow-scale Corg distribution. Strong correlations between Corg, grain size, and edge proximity suggest that current attenuation increases fine-sediment deposition and, therefore, carbon burial with distance into the meadow. By mapping the sediment Corg pool, we provide the first accurate quantification of an enhanced carbon stock attributable to seagrass restoration. The top 12 cm of the bed contain 3660 t Corg, approximately 1200 t more Corg than an equal area of bare sediment. Most of that net increase is concentrated in a meadow area with low tidal current velocities. Managers should account for the effects of meadow configuration and current velocity when estimating seagrass blue carbon stocks. Our results suggest that a large, contiguous meadow should store more blue carbon than an equal area of small meadow patches.


Introduction
Seagrass meadows are highly productive ecosystems [1] that bury organic carbon [2,3], making them sinks in the global carbon cycle [4]. Bed accretion from canopy particle trapping, seston PLOS  shear stress to a constant level over short distances (perhaps <0.5 m, cf. [33,34]), the 'edge zone' may be narrow and have a negligible impact on the total C org stock ( Fig 1B). Alternatively, external sediment supply may decrease with distance as particles settle out [32], causing sediment accumulation to peak and then decrease towards the meadow center ( Fig 1C). Canopy filtration of external particulate organic matter [35] and deflection of incoming wave energy [36,37] may result in more C org burial near the edge. Some empirical evidence now exists for seagrass sediment C org concentration variability at patch-and regional-scales [24,27], but there is a lack of empirical support to link site and process-based studies with spatial patterns at the meadow-scale. Few seagrass studies sample over whole meadow areas and most do not account for spatial autocorrelation (for example, [24]), thereby limiting their ability to identify statistically significant controls on C org accumulation at spatial scales !1 km.
The existence of meadow-scale C org spatial gradients would have implications for blue carbon finance for seagrass restoration. A framework now exists for allocating carbon dioxide equivalent (CO 2 e) offset-credits to restoration projects. This framework, Verified Carbon Standard methodology VM0033 [38], requires that managers stratify project areas to account for spatial variability when estimating blue carbon stocks (cf. [39] for information on stratification). If zones of higher and lower C org accumulation exist within meadows, managers must identify these zones (i.e. strata) and scale average C org concentrations within each zone separately to generate a whole-meadow stock estimate.
The restored Zostera marina (eelgrass) meadow in South Bay, Virginia, provides a unique opportunity to investigate sediment C org spatial patterns at the meadow-scale. This mature meadow has a well-documented restoration and expansion history [40,41,42,43], allowing us to consider C org accumulation time (i.e. site meadow age) as an independent variable potentially affecting the C org spatial distribution. Eelgrass seeds were broadcast in South Bay in 0.5 and 1-acre seed plots from 2001-2004 [40], which later coalesced into a single, contiguous meadow encompassing >6 km 2 [43], making it the world's largest restored seagrass meadow. This meadow is accumulating C org from both eelgrass and allochthonous sources relative to bare sites [3,42,44]. Wind-generated shear stress is a dominant control on sediment suspension throughout this coastal bay system [45,46]. The meadow attenuates both wave and tidal currents, thereby affecting suspended sediment concentrations and sediment accretion rates relative to bare areas [3,47]. We expected to find sediment C org spatial gradients related to edge proximity, in addition to a C org distribution pattern related to meadow expansion. Spatial gradients related to differences in relative location and meadow age may be attributable to covarying factors known to affect sediment accumulation at individual meadow sites, especially shoot density [26,42].

Study area
The restored South Bay eelgrass meadow and surrounding area is part of the Virginia Coast Reserve, a Long-Term Ecological Research site (VCR-LTER). South Bay is separated from the Atlantic Ocean by Wreck Island to the east and bordered to the west by Man and Boy Channel. The meadow bathymetry is shallow, with an average depth at MSL of 0.76 m and a standard deviation of 0.28 m [48]. Tides enter and exit the area via inlets to the north and south of the meadow. The mean tidal range is 1.2 m (Fig 2) [46]. The Virginia coastal bays are oligotrophic [49], resulting in a light environment conducive for plant growth [42,50]. Shoot densities exceed 600 shoots m -2 at maximum biomass in the summer, and the canopy height ranges from 32-45 cm [42]. The east half of the meadow adjoins Wreck Island; the west half of the meadow has a well-defined edge that separates the meadow from bare, subtidal areas (Fig 2). South Bay's large surface area, relatively constant bathymetry, and a water depth to canopy height ratio of approximately 2:1 at MSL make it an ideal location for assessing canopy-hydrodynamic effects. Restoration has resulted in increased sediment carbon concentrations in the meadow [42], such that a mid-meadow site exhibited significantly higher carbon content than an adjacent bare site in 2012, 0.52±0.010% versus 0.36±0.012% [3]. Hansen and Reidenbach [30,47] measured bed shear stress and suspended sediment concentrations at locations both inside and outside of this meadow from 2010 to 2011. Average Reynolds stresses were lowest in the southwest part of the meadow, approximately 0.5 cm 2 s -2 (z = 0.5 m), compared with~1 cm 2 s -2 at a mid-meadow site and 1.5 cm 2 s -2 outside of the meadow [30]. Suspended sediment concentrations averaged~20 mg l -1 at a mid-meadow site and~30 mg l -1 outside of the meadow [47].

Data collection and sample preparation
Sediment samples were collected along eight subtidal transects in July 2013 to quantify C org distribution patterns in the meadow (Fig 2). Parallel transects were laid out from meadow edge to interior in each cardinal direction to provide broad meadow coverage. The transect sites were arrayed systematically (note that permits were not required for sediment collection from public bottomlands). Each transect contained eight sites spaced 150 m apart with the exception of transect 1, which had ten sites. Eight sites fell within original restoration seed plots [40] and the others represented meadow ages ranging from <1 to 12 years due to natural meadow expansion. Four sites near the meadow provided a bare control group (Fig 2). Four replicate 60 cc hand cores were collected at each site to a depth of 12 cm and divided into four 3-cm intervals to generate a sediment C org profile for each site. The bed has aggraded~3-4 cm due to restoration [3]. Each of the 264 cores collected during this study were visually inspected for compaction when collected, which was approximately 7%, given the predominantly sandy sediment in this study area (mean grain size = 71 μm [45]).
Macroscopic roots, rhizomes (i.e. belowground biomass), and shell fragments were removed from samples to isolate the sediment organic matter (OM) component from belowground biomass. Sediment bulk density, %OM, and percent carbon and nitrogen (%C and % N) were determined following standard methods used previously at this site [3,42]. Loss on ignition (LOI) in a muffle furnace at 500˚C for six hours was used to determine %OM. A Carlo Erba NA 2500 Elemental Analyzer was used to determine %C and %N. Bulk C org was determined using the element analyzer method described in [20].
Meadow age at each site was established using aerial photographs taken annually beginning in 2001 (high resolution images provided by R. Orth, VIMS [51]). Sample sites were georeferenced relative to each aerial photo in ArcGIS 10.2. By observing meadow presence/absence in each photo, we determined the length of time that seagrass has been present with one-year precision. The original seed plots coalesced over time into a single meadow, which continued to expand naturally, such that seagrass remained present at every meadow site after its first appearance (Fig 2).
The georeferenced aerial photographs also allowed us to determine site distance from the meadow perimeter. Site edge proximity was measured two different ways: linear distance along transects to the 2013 meadow perimeter and Euclidean distance to the 2013 edge (Near analysis in ArcGIS 10.2). The first measure allowed us to compare C org concentration changes with distance along a given transect. The second measure established site location relative to the meadow boundary that intercepts incoming current and wave energy.
Several additional variables were also measured that could influence the C org spatial distribution pattern, including site relative exposure, peak-summer (July) shoot density (shoots m -2 ), sediment grain size distribution (top 3 cm of the bed), mean water depth (bathymetry), and sediment C:N ratio. A relative exposure index (REI) was calculated for every site according to methods in [21]. Effective fetch was found by intersecting radiating lines at each site with surrounding land surfaces delineated using aerial imagery in ArcGIS 10.2. Wind vector and frequency data necessary for REI was obtained from a LTER monitoring station immediately south of the meadow on Godwin Island [52]. Replicate shoot density counts and grain size samples were collected at a subset of sites that were randomly selected to provide broad spatial coverage (n = 16). Average density counts were also taken at six additional, mid-meadow sites during a VCR-LTER annual survey. Shell fragments were removed from grain size samples, which were oxidized using 30% H 2 0 2 and then acidified using a 5.0 pH acetic acid and NaOAc solution to remove OM and carbonates. Grain size distributions were determined using a Beckman-Coulter BLS 13 320 Laser Diffraction Particle Size Analyzer. Water depth at each site was determined by extracting bathymetric data [48] by site in ArcGIS 10.2. REI and bathymetry provided measures of relative wave energy and tidal current strength, respectively [53]. Mean grain size provided a measure of time-integrated shear stress and water residence time [54,55]. The C:N ratio provided an indication of OM source.

Analyses
We identified meadow-scale spatial patterns by kriging sediment %OM, bulk C org , bulk C:N, shoot density, and grain size (mean grain size and % sand fraction) in ArcGIS 10.2 Geostatistical Analyst. Kriging accounted for spatial autocorrelation and provided error estimates for interpolated values (cf. [56] and references therein), which allowed us to generate robust distribution maps. We fit multiple variogram models to each dataset (stable, circular, spherical, exponential, and Gaussian) and cross-validated the models using root mean square errors. Summing the kriged C org distributions allowed us to quantify the gross C org stock to a depth of 12 cm. We subtracted the average background concentration from each map cell by bed interval to assess the enhanced (i.e. net) stock attributable to the meadow.
We evaluated observed gradient relationships between C org , meadow age, and edge proximity using regression analysis (lm, stats package) and Kendall correlation (rcor.test, ltm package) in R [57] to determine whether the variables met multiple regression assumptions. Correlation analysis allowed us to check for possible multi-collinearity between age and edge proximity (Euclidean distance). Linear regression analyses allowed us to determine whether variable relationships exhibited linearity-a multiple regression assumption. Shapiro-Wilks tests were run on C org concentrations by depth interval (shapiro.test, stats package) to verify that the response variables were normally distributed. In addition to age and edge proximity, our primary variables of interest, we also considered gradient relationships between C org and several potential explanatory variables: REI, density, mean grain size, and bathymetry. The Kendall correlation analysis included kriged values for density and grain size, in addition to measured values, to obtain equal sample sizes for all six independent variables. Spatial autocorrelation potentially confounds attempts to determine the relative importance of age and edge proximity in predicting sediment C org concentrations at the meadow-scale. To address this, we compared meadow age and edge proximity (Euclidean distance) effects on sediment C org using spatial autoregressive models (spdep package version 0.6-4 for [57]), which utilized a neighborhood weights matrix and Moran's I to account for autocorrelation (for example, [58]). We evaluated the C org data for both the 0-3 cm and 3-6 cm depth intervals. The top~3 cm of the bed accumulated after the start of the restoration [3]. The 3-6 cm interval approximately corresponded with the rhizosphere in this system. Spatial lag and error dependencies were evaluated using Lagrange Multiplier tests (lm.LMtests, spdep package, see [59]). The best models were determined by maximum likelihood estimation using Akaike's Information Criterion (AIC).

Spatial distribution of sediment C org
Sediment C org concentrations varied across the meadow and with sediment depth. The meadow-wide gross concentration for the top three cm of the bed averaged 3.92±0.15 (SE) mg C org cm -3 ; the 3-6 cm bed depth interval had the highest average concentration, 5.66±0.25 (SE) mg C org cm -3 ( Table 1). The gross concentration ranged from 1.42 to 7.19 mg C org cm -3 in the top three cm and from 1.66 to 9.84 mg C org cm -3 in the 3-6 cm interval. However, this variability was non-randomly distributed across the meadow. C org concentrations showed significant spatial autocorrelation at all distances 1160 m for all bed depth intervals (Moran's I > 0.04, p < 0.05), resulting in strong spatial gradients. Sediment C org decreased along each of the six transects extending from the meadow interior to the edge (Fig 3 and Table 2). In contrast, C org concentrations increased with distance from the meadow interior to the perimeter on the two transects adjacent to the barrier island (T5 and T7 in Fig 3).
Interpolating the %OM, bulk C org , and C:N data consistently yielded two discrete spatial regimes: a kriged zone of higher %OM and C org encompassing most of the southeastern half of the meadow and another zone of decreasing %OM and C org to the northwest (Fig 4). C:N ratios and grain size data yielded a similar kriged pattern, with higher C:N ratios and a higher percentage of larger grains to the northwest and lower values to the southeast, near Wreck Island (Figs 4 and 5). Summing the interpolated sediment C org within the meadow area to a bed depth of 12 cm gave a total meadow (gross) stock of 3662 t C org . Subtracting average background concentrations measured at the bare sites from each meadow site by depth interval and interpolating the net increase gave a net stock of 1173 t C org .

Spatial variables
Several factors possibly account for the spatial distribution of sediment C org . Edge proximity, meadow age, grain size and bathymetry were all significantly correlated with C org at the meadow-scale. Edge proximity regression relationships (Euclidean distance) were highly significant at all four depth intervals. The highest adjusted-r 2 values were for the 0-3 and 3-6 cm intervals ( Table 2). C org and meadow age were highly correlated, but exhibited a relatively weak, positive linear regression relationship (Tables 2 and 3). The strongest regression and correlation relationships were between the 0-3 cm C org and sediment grain size (Tables 2 and  3). C org concentrations were not significantly correlated with shoot density or REI (Table 3). Unlike the kriged grain size distributions, the kriged density distribution did not match the C org distribution (Fig 5).
Edge proximity, age, density, grain size, and bathymetry were also significantly correlated with one another, potentially indicative of landscape-scale interactions among the variables. Grain size was negatively correlated with edge proximity, age, and bathymetry, and positively correlated with REI ( Table 3). The strongest correlation coefficient was between grain size and site distance from the open perimeter (τ = -0.575). Site age and edge proximity (Euclidean distance) were also moderately correlated but not co-variates (τ = 0.574), allowing us to compare  Seagrass blue carbon spatial patterns their ability to predict C org concentrations using multiple regression. The meadow has expanded outward over time, but it also coalesced in places and filled in behind Wreck Island relatively recently (Fig 2). Several factors complicated our ability to compare the effects of edge proximity and age on C org concentrations directly through multiple regression. C org varied non-linearly with age (Fig 6), so we applied a log-log transformation, which yielded linear relationships with adj-r 2 values >0.3 for both the 0-3 and 3-6 cm data ( Table 2). The log-transformed C org data for the 0-3 cm interval met the dependent variable normality assumption (W = 0.972, p = 0.143) and also varied linearly with the edge data. The log-transformed C org data for the 3-6 cm interval was not normally distributed (W = 0.963, p = 0.044). Running a multiple regression analysis comparing edge and log-transformed age on C org concentrations within the top three cm of the bed yielded spatially autocorrelated regression residuals (Moran's I = 0.535, p<1.10E-4). Lagrange Multiplier tests found strong spatial lag dependence (RLMlag = 6.1796, p<0.013). We accounted for this spatial autocorrelation using spatial lag and error analyses, both of which identified edge proximity-not age-as a significant predictor of C org concentrations at the meadow-scale (Table 4).

Meadow-scale controls on sediment C org accumulation
Our results indicate that differences in relative location within the meadow affect the C org stock distribution and overshadow other factors, including seagrass age and shoot density, that are known to affect C org concentrations at the plot-scale. We show that edge proximity affects C org concentrations over much larger spatial scales than previous studies have recognized, potentially resulting in seagrass meadow spatial gradients >1 km in length. Carbon stock estimates should take these potential meadow-scale spatial patterns into account.  Rather than reflecting the meadow's expansion history over the preceding 12-year period, the meadow-wide C org distribution appears broadly consistent with the hypothesis that current attenuation promotes higher C org concentrations with distance from the edge (Fig 1A). The meadow-wide grain size distribution (Fig 4) also supports this hypothesis. Previous studies have noted that suspended sediment is advected into this meadow [30], increasing the percentage of silts and fine sands at meadow sites over time [42]. This suspended sediment fractionates according to particle size as it is deposited across the meadow, with finer particles settling out in the southeastern meadow, where Hansen and Reidenbach [30] documented the lowest average Reynolds stress. This, in turn, facilitates more C org storage, because smaller grains have more surface area for C org adsorption [55]. This process appears to be driven largely by tidal currents. Hansen and Reidenbach [30] observed similar Reynolds stresses attributable to wave-dominated flows in different meadow locations, which may explain why we did not observe a significant correlation between C org and REI. In addition to canopy current attenuation, basin geomorphology possibly accounts for some of the reduction in current velocity in the area of the meadow adjacent to Wreck Island and furthest from the two inlets. We note that a C org spatial pattern is weakly present within the underlying 9-12 cm depth interval (Fig 4), which was deposited prior to the meadow restoration [3]. However, the meadow has     Table 2 for individual regression statistics.
clearly accentuated the observed pattern, as evidenced by the magnitude of the discrepancy in C org concentrations between the two spatial regimes within the top 6 cm of the bed (Fig 4). Root and rhizome-derived carbon compounds may also contribute to sediment C org accumulation below 6 cm. These results confirm that large, sediment C org spatial gradients are possible and should be considered when estimating blue carbon stocks. Similar studies are now needed to determine how varying current velocity, meadow configuration, and the water depth to canopy height ratio might give rise to particular sediment C org gradients at this spatial scale. Other C org distributions may be possible in other meadows. For example, an even larger meadow with a lower water depth to canopy height ratio might experience lower C org concentrations further from the edge, because of reduced sediment delivery ( Fig 1C). However, given that we observed increasing sediment C org concentration gradients (Fig 1A) >1 km in length, and given that smaller, patchy meadows are relatively common, edge proximity possibly limits C org accumulation in many systems. Other studies have speculated along these lines but lacked the ability to control for meadow expansion as a possible confounding variable (for example, [27]). Our results also highlight the importance of considering spatial autocorrelation and its potential effect on measured quantities at individual sites within a given seagrass meadow.
The potential importance of external sediment raises the possibility that much of the sediment C org stored in this meadow is, in fact, allochthonous in origin. Greiner et al. [44] found that only half of the sediment C org at an interior site in the South Bay meadow derived from vascular plants. The C:N ratio in the C org hotspot in the southeast part of the meadow (Fig 4) conforms more closely to the range observed for seston and macroalgae than for Z. marina in this system [60]. The high C:N values in the northwest half of the meadow more closely resemble Z. marina. The grain size spatial distribution suggests that the southeastern meadow experiences lower current velocities and longer residence times, both of which possibly increase seston accumulation, which would increase the magnitude of observed C org spatial gradients across the meadow. However, additional isotopic work is needed to conclusively identify C org sources at this spatial scale. The fact that edge proximity, not age, significantly predicts meadow-wide sediment C org concentrations indicates that meadow-and regional-scale factors should be considered when estimating whole-meadow carbon stocks. Recent studies consider blue carbon accumulation as a function of plot-scale factors, including meadow age and plant density (for example, [61]), without considering possible spatial scale effects. Age and C org concentrations are positively correlated at individual sites in this study, but differences attributable to relative location overshadow differences due to age at this spatial scale. Shoot density also affects sediment C org accumulation at the plot-scale [42,62], but density alone is not a good proxy for site C org concentrations at the meadow-scale. This is likely due to the fact that meadow canopy structure (i.e. shoot density and biomass) varies considerably over small spatial and short temporal scales [42,63] and because density effects on sediment resuspension appear to be non-linear [45]. Consequently, a snapshot assessment of canopy structure would not necessarily correspond with the sediment C org distribution, which reflects the balance of accumulation and resuspension over interannual timescales. Likewise, REI might correlate with C org in isolated seagrass patches but does not account for current attenuation by the canopy.

Implications for financing seagrass restoration using blue carbon offsetcredits
The distribution of sediment C org in this meadow follows approximately linear C org concentration gradients-not irregular zones of higher and lower C org concentrations controlled by wind fetch, canopy complexity, age, or other factors. A representative C org concentration for stock estimation might, therefore, be obtained by averaging samples collected from a relatively small number of sites distributed along the gradient. However, managers should avoid overestimating C org stocks by relying on point-based literature values, models, or default values-all permissible approaches under VM0033 [38]. Scaling the C org measurement reported for this meadow by Greiner et al. [3] by the total meadow area would overestimate the gross sediment C org stock by almost 20%, because of the spatial gradients. Likewise, managers should not simply scale C org accumulation model results calculated for small unit areas (for example, [61]) or rely exclusively on near-surface 210 Pb accumulation rates that do not account for remineralization [64], without understanding possible meadow-scale effects.
Regarding carbon offset-credit finance for seagrass restoration, the enhanced sediment C org stock attributable to the meadow after more than a decade translates to approximately 4,300 t CO 2 . Incorporating sequestered Z. marina biomass C org would increase this total [12]. However, the aboveground biomass is sloughed off, and the fate of this exported C org is uncertain, so the sequestered stock would correspond to the annual cycle average, not peak standing biomass.
Managers might be able to increase blue carbon storage by considering meadow configuration, basin geomorphology, and regional hydrodynamics when locating seagrass restoration sites. A large, contiguous meadow should store more blue carbon than an equal area of small meadow patches. If blue carbon storage is a management goal, restoration should be initiated at sites that are suitable for the accumulation of fine sediment. As the meadow expands, these locations should accumulate more blue carbon, due to scale-dependencies observed in this study, and adjacent areas should begin to bury blue carbon, thanks to the canopy particle-trapping feedback [29,37]. Additional studies that determine how current velocity, meadow configuration, and water depth interact to influence meadow-scale C org gradients can aid blue carbon accumulation modeling efforts at spatial scales relevant to restoration managers.

Conclusions
This study indicates that edge proximity can better explain a seagrass meadow's sediment C org distribution than spatial differences in accumulation time. Although meadow age and seagrass shoot density affect C org accumulation at the plot-scale, these drivers can be overshadowed by differences in relative location at the meadow-scale. Progressive canopy attenuation of currents may explain the C org distribution observed in this study. As currents move through the canopy, suspended sediment becomes stratified and is deposited according to particle size, which likely facilitates more C org burial at more interior sites, irrespective of site meadow age. These findings highlight the potential importance of external sediment for seagrass blue carbon accumulation and the need to consider meadow-scale spatial gradients when quantifying whole-meadow carbon stocks.