Sampling took place in 20 seagrass meadows in the western part of the Baltic Sea, along the coasts of Schleswig-Holstein (n = 17) and Mecklenburg-Vorpommern (n = 3) in northern Germany ( Figure 1 ). The German Baltic Sea coast consists of shallow bays and fjords that experience weak water currents and low wave heights (Petterson et al., 2018). Like the whole Baltic Sea region, German coastal waters were geologically shaped by the last glacial periods (Schmölcke et al., 2006; Andrén, 2012). As a consequence of the last deglaciation, a conglomerate of differently sized stones, sand, and clay settled in the southwestern part of the Baltic Sea basin. The seabed consists of shallow sandy and muddy layers with consolidated marl underneath. A maximum water depth of 40 m is reached in the western part of the Baltic Sea. However, seagrasses here are rarely observed deeper than 8 m seawater depth (Schubert et al., 2015). The Baltic Sea is the largest brackish water basin in the world and because of the narrow Danish Straits connecting the Baltic Sea to the North Sea, low rates of water exchange are observed (residence time of 35-40 years) resulting in high eutrophication from nutrient discharge by the nine Baltic Sea nations that enclose it (HELCOM, 2018). However, eutrophication, which negatively impacts seagrass health and bathymetric range, is most pronounced in Germany, Russia and Poland (Thorsøe et al., 2022).

Sites in the study area represent the greatest environmental gradient, ranging from wave exposed (e.g. Heidkate, Falshoeft lighthouse, Teichhof, Goehren; Figure 1 star shape) to sheltered (e.g. Orth, Maasholm; Figure 1 square shape), relatively pristine (nature reserve e.g. Gelting Bay and Graswarder; Aschau is protected from human impact due to military controlled access to the area; Figure 1 diamond shape) to varying degrees of anthropogenic inputs ( Figure 1 circle shape), such as adjacent to a marina (e.g. Gluecksburg, Glowe; Figure 1 pink circles); heavy ship traffic and near urban areas (e.g. Kiel Fjord, including Falckenstein, Seebar, Hasselfelde; Figure 1 purple circles); agricultural land (e.g. Wackerballig; Figure 1 orange circles); tourist areas (e.g. Grossenbrode, Kellenhusen, Sierksdorf; Figure 1 red circles), and proximity to river influx (e.g. Gahlkow, Hasselfelde, Niendorf; Figure 1 blue circles).

A total of 169 sediment cores (30 cm length; 5.5 cm inner diameter) were sampled in seagrass meadows (n = 110 cores) and nearby unvegetated sediments (n = 59 cores). Nine cores (with some exceptions) were collected at each site, from three sublocations: (1) in the high-density part of the meadow (‘dense seagrass’ hereafter), (2) low density or fringe of the meadow (‘sparse seagrass’ hereafter), and (3) adjacent unvegetated sediments at least 5 m from seagrass (‘unvegetated’). Dense and sparse seagrass sublocations are collectively referred to as ‘seagrass-vegetated’ sediments or sublocations. Sediment cores taken within the seagrass meadow were sampled at least 10 m apart from each other. Seawater depth was measured from a diver computer. Seagrass shoot density was counted using a 0.04 m² frame, and the leaf lengths of five randomly selected plants were measured next to the location where the core was extracted. Seagrass density and leaf length vary between sites and at meadow stage of maturity. The measure of “seagrass complexity” was defined as the product of seagrass canopy height and shoot density, to obtain the sum of leaf heights within a unit of area (in m/m²) (Prentice et al., 2019). Seagrass canopy height and shoot density exhibit seasonal variations, and so it must be noted that all sampling took place between July 15^th and September 1^st of 2019 or 2020, at the peak of the seagrass growing season for Z. marina in Germany, except for the samples in Gelting Bay and Falshoeft lighthouse which were collected in late November (23^rd and 24^th of 2020) due to permitting and accessibility issues.

Cores were collected between 1-5 m seawater depth manually via SCUBA divers (self-contained underwater breathing apparatus) pounding impact resistant PVC tubes into the sediment with a rubber mallet. Cores were capped at both ends by the divers and stored upright for transport to shore, and in a cooler thereafter for transport to the lab. Cores were stored at 0°C until further processing.

Sediments are known to shift during the coring process, so the degree of compression was measured by divers once the core was fully inserted. Compaction was calculated as the distance (in cm) from the top of the core to the sediment surface outside of the core, divided by the sample depth (cm of compression/cm core depth). A compression correction factor was calculated by dividing the length of the sample recovered by the length of core penetration (Howard et al., 2014). Cores were subsectioned into 5 cm depth intervals, and depth intervals were adjusted based on compaction.

Sediments from each depth interval were homogenized, measured for dry-bulk density and subsampled for chemical analyses (total C_org, δ¹³C, δ¹⁵N, ¹⁴C, each described in detail in subsections below). A total of 394 sections underwent total C_org analyses: all top (0-5 cm) layers of all cores were processed for total C_org, but the selection of subsequent core depth intervals depended on color changes throughout the core (color changes indicate potential changes in C_org content) (Nederbragt et al., 2006, and e.g. Klingler et al., 2020). If no color change was apparent, the mid (10-15 cm) and bottom sections of the core were processed (typically 20-25 cm, except where cores were shorter due to marl or difficult sediments, such as large rocks or bolder field beneath the surface of the sediment). If a color change occurred, the respective section was instead examined. It must be noted that the total C_org was measured in 0-5 cm, 10-15 cm and 15-20 cm sections of all sites, but not all cores as described above. The 5-10 cm depth interval of all cores was used for grain size analysis (described below, under ‘grain size analysis’). Visible plant material (seagrass roots and rhizomes) and infauna (lug worm Arenicola marina, and, in some cases, soft tissue of clams and mussels) were removed from the sediment prior to chemical analyses (hereafter referred to as ‘visible organic’), but all other materials (e.g. such as wood) were left in the sediment (referred to as ‘sediment organic carbon’ or ‘SOC’). Shells, were also left in the sediment after the organic tissue was removed.

These values were represented as maximum wave-generated orbital velocity (MOV) at the seafloor, modelled for the year 2021. Due to the limited availability of simulated wave data, MOVs were calculated only for cores taken at the outer coast sites of Schleswig-Holstein, not for the Schlei fjord site (Maasholm) and also not for any of the sites in Mecklenburg-Vorpommern (Gahlkow, Glowe, Goehren). MOV was calculated as a function of wave height, mean wave period, mean wave length (all simulated values), and seawater depth for each site according to linear wave theory as per the procedure described in Bobsien et al. (2021). MOV averages (‘avg. MOV’) and maximum MOV (‘max MOV’) were based on MOV calculations for simulated hourly seastate values within the year 2021 (n = 365 days x 24 hours). Average maximum MOV (‘avg. max MOV’) at each site was obtained by averaging the highest MOV observed at the depth and location of each core.

The 5-10 cm depth interval of each core was processed for grain size distribution. Visible organics were first removed, then the sediments were oven dried at 60°C for a minimum of 48 hrs. Dry sieving was conducted with stainless steel Test Sieve ISO 3310-1 and sieve shaker (Fritsch Analysette 3 Spartan Pulverisette 0) set to amplitude 2 mm for 15 min. The fractions of sediment in 2 mm, 1 mm, 500 μm, 250 μm, 125 μm, 63 μm, and <63 μm (herein referred to as ‘fine sediment’ fraction) size classes were determined to the nearest 0.01 g, and the percent amount of each size class was calculated using the total sample mass obtained from the sum of each fraction.

Sediment total C_org was determined using an Elemental Analyzer (EURO EA Elemental Analyzer). Sediments were first dried at 60°C for 48 h and then ground to a homogeneous fine powder using a mechanical agate ball mill (Fritsch Pulveisette 5) set at rotational speed 240-300 rpm for 15-20 minutes. A subsample of this homogenized sediment was acidified to remove inorganic carbon by adding 1 M HCl drop-by-drop until gas evolution ceased. Samples were observed under a dissecting microscope to ensure CO₂ had fully evolved. These acidified samples were re-dried at 60°C for 48 h, ground again (with mortar and pestle), and encapsulated into silver capsules. Total organic carbon concentrations were calculated based on a linear regression. Acetanilide and sediments were used as standards to measure data accuracy and analytical uncertainty. For sections that were not measured by EA, the section carbon content was taken from the quantify produced in the layer above (with similar coloration, see core subsampling method above for further details).

To decipher the source material of the remaining (non-visible) organic fraction in the sediments (referred to as SOC), four sites (Falckenstein, Graswarder, Grossenbrode, Wackerballig) underwent further examination using biotracers of δ¹³C and δ¹⁵N. Samples were pre-treated as outlined above (see ‘C_org stock quantification’), and then combusted in an elemental analyzer system (NA 1110, Thermo) coupled to a temperature-controlled gas chromatography oven (SRI 9300, SRI Instruments), connected to an isotope ratio mass spectrometer (DeltaPlus Advantage, Thermo Fisher Scientific) as described by Hansen et al. (2009). δ¹⁵N, δ¹³C ratios were reported in delta notation relative to the international Air-N₂ and VPDB scale following the equation:

where R = ¹⁵N/¹⁴N or ¹³C/¹²C. N₂ and CO₂ gases were used as reference gases and calibrated against International Atomic Energy Agency (IAEA) reference standards (N1-,N2-, NO3-) and National Institute of Standards and Technology (NBS-22 and NBS-600) compounds. To measure analytical uncertainty, acetanilide and caffeine were used as internal standards after every sixth sample to test if the analytical setup was working properly. Precision for Caffein was ± 0.13 ‰ for δ¹⁵N and ± 0.09 ‰ for δ¹³C; and ± 0.20 ‰ for δ¹⁵N and ± 0.27 ‰ for δ¹³C for Acetanilide.

A two biotracer (δ¹³C and δ¹⁵N), five-source Bayesian mixing model using R package (R Core Team, 2022) MixSIAR (Stock and Semmens, 2016) was used to determine the contribution of C_org to seagrass-vegetated and unvegetated sediments in Falckenstein, Grossenbrode, Graswarder, and Wackerballig. Signatures of five potential sources were compiled from previous studies in Kiel Fjord, Germany and Gulf of Gdansk, Poland, and consisted of: (1) Pilayella littoralis, a filamentous brown algae known to form thick drifting mats accumulating in seagrass meadows in Germany (Kruk-Dowgiallo, 1991; Kiirikki and Lehvo, 1997), (2) other macroalgae (e.g. Ulva intestinalis and Cladophora fracta) extracted from Table 2 in Maksymowska et al., 2000; and (3) phytoplankton (collected by seston), (4) epiphytes attached to seagrass leaves, (5) seagrass leaves (there is little difference between above- and below-ground biomass signatures, see Röhr et al., 2016) extracted from Table 1 in Mittermayr et al. (2014). Sampling location (four locations) and vegetation state (two states: seagrass-vegetated vs unvegetated) were included as random effects in the stable isotope mixing model. Vegetation state was nested within sampling location for the random structure. We specified an uninformative, generalist prior and discrimination factor of 0 as no further fractionation is expected in stored C_org (Greiner et al., 2016; Miyajima et al., 2017). Note that visible seagrass roots and rhizomes were omitted from this analysis, which may lead to a lower contribution of autochthonous organic carbon sources.

Large amounts of exceptionally well-preserved wood pieces were discovered in sediment cores extracted from the seagrass meadow in Sierksdorf, Luebeck Bay ( Figure 2 ). Of those, eleven pieces from a depth of 5 to 24 cm in four sediment cores were radiocarbon dated. The samples were visually inspected under a microscope and an appropriate amount of wood material was selected for dating. A standard decontamination procedure was subsequently applied to remove carbonates and soil humic contaminants, consisting of 1% HCl, 1% NaOH at 60°C and again 1% HCl. All radiocarbon measurements were conducted at the Leibniz-Labor, using the type HVE 3MV Tandetron 4130 accelerator mass spectrometer (AMS). Following standard procedures, the ¹⁴C/¹²C and ¹³C/¹²C isotope ratios were simultaneously measured by AMS, compared to the CO₂ measurement standards (oxalic acid II), and corrected for effects of exposure to foreign carbon during the sample pretreatment. The resulting ¹⁴C-content was corrected for isotope fractionation, related to the hypothetical atmospheric value of 1950, and reported in pMC (percent Modern Carbon). This value was used to calculate the radiocarbon age according to Stuiver and Polach (1977). The reported uncertainty of the ¹⁴C result takes into account the uncertainty of the measured ¹⁴C/¹²C ratios of sample and measurement standard, as well as the uncertainty of the fractionation correction and the uncertainty of the applied blank correction. The radiocarbon ages were translated to calendar ages using the software package OxCal4 (Ramsey and Lee, 2013) and the Intcal20 dataset (Reimer et al., 2020). Note that this is not an observation unique to Sierksdorf as our diver and drop camera surveys showcase old wood material that has become visible on the fringes of seagrass meadows where erosion has occurred: in other parts of Sierksdorf that were not cored in this study, Neustadt, Fehmarnsund, Kellenhusen, Heidkate, wendtorf, Dollerupholz ( Table 1 ).

Statistical analyses were performed in R version 4.1.3 (R Core Team, 2022) for Mac OS X, and all mixed models were implemented with package “lme4” (Bates et al., 2015). For all tests, significance was determined at p < 0.05. Visual inspection of standard model validation graphs was used to verify model assumptions: residuals versus fitted values were used to verify homogeneity; a histogram or Quantile–Quantile (q–q) plot of the residuals for normality; and residuals versus each explanatory variable to check independence.

Generalized linear mixed models (GLMM) with a Gamma distribution and log link function were used to test: (1) core position effect on C_org between sublocations (dense seagrass, sparse seagrass, unvegetated) within each site (locally), with a random intercept of sublocation nested in sampling location (site), and (2) biophysical factors influencing regional differences in stocks within seagrass-vegetated sediments, with sampling location as a random intercept. In the local model (1), a posthoc test using Tukey contrasts was performed with package “multcomp” to examine multiple comparisons between sublocations of each site. A multi-model inference approach based on AICc was used to determine the strongest predictors of regional C_org (2). Selection criteria AICc was used (instead of AIC) due to a small sample size (n/k < 40; n = total number of observations; k = total number of parameters in the most saturated model, including both fixed and random effects). The saturated model included four factors: seagrass complexity, average MOV, seawater depth, percent fine sediments, and their interactions. However, it is important to consider multicollinearity when interpreting model outputs as coefficient estimates (i.e. beta coefficients) and p-values become very sensitive to any small changes in the model. Multicollinearity of parameters in the dataset and among saturated model terms were examined using the Pearson correlation coefficient (ρ) and variance inflation factor (VIF). Correlation coefficients < 0.6 and VIFs approx. < 5 indicate an acceptable level of collinearity, VIFs > 10 warrant further investigation (Montgomery and Peck, 1992). Low collinearity was observed in the data itself (ρ = -0.43 to 0. 35), but structural multicollinearity (between some model terms) was observed in the saturated model (VIF 1.876 to 14.051), so all terms with VIF >10 were removed from the global model. Hence, the final, global model included four factors: seagrass complexity, average MOV, seawater depth, percent fine sediments and two-way interactions between all variables except average MOV and fine sediment fraction, and all three- and four-way interactions were excluded (excluded VIFs: 10.055 to 14.051). After model selection, statistically indistinguishable (ΔAICc <2) candidate models were averaged using package ‘MuMIn’ (Bartoń, 2022). To estimate the importance of each variable, the RVI (relative variable importance) value was computed by summing Akaike weights across all averaged models where a particular variable appeared. Model-averaged coefficients and p-values were obtained from the ‘full average’ (rather than the ‘conditional average’).

Averaged across all sediment cores in seagrass-vegetated sublocations, including dense and sparse seagrass sublocations (n = 110 cores total) and integrated to 25 cm core length, C_org stocks averaged at 1,920 ± 402 g C/m², and varied 22-fold between sites, ranging from 475 ± 61 and 10,577 ± 6,178 g C/m² ( Figure 3 ). For seagrass-vegetated sediments, the highest average C_org stocks were found in Maasholm (Schlei Fjord) and Sierksdorf (Luebeck Bay), while the lowest values were observed in Heidkate (Kiel Fjord) and Grossenbrode (Luebeck Bay).

An examination of the local (within site) variation in C_org content across all sites showed significantly greater (by three-fold) stocks in seagrass-vegetated sediments (both dense and sparse seagrass) compared to nearby unvegetated sediments (dense seagrass: average 1,523 ± 244 g C/m²; sparse seagrass: 2,316 ± 795 g C/m²; unvegetated: 611 ± 93 g C/m²) ( Table 2 ; Figure 3 ). However, considering sites individually, within site comparisons ranged widely: in all but four cases (Gluecksburg, Gohren, Hasselfelde, Seebar), seagrass-vegetated sediments had 1.5 to 31 times more C_org than nearby unvegetated sediments, with Sierksdorf and Maasholm showing the greatest discrepancies (10 and 31 times that of unvegetated sediments). Gluecksburg and Hasselfelde had the same; Gohren, and Seebar less C_org than nearby unvegetated sediments. Seebar had the highest average C_org stocks, twice as much as in nearby seagrass-vegetated areas. C_org content did not differ significantly between sparse and dense seagrass sublocations.

Overall, the stable isotope signatures of the SOC fraction of seagrass-vegetated sediments (including dense and sparse seagrass sublocations along the entire core depth) were near identical to those of unvegetated sediments (δ¹⁵N 5.2 ± 0.3‰, δ¹³C -21.9 ± 0.2‰ vs δ¹⁵N 5.3 ± 0.2‰, δ¹³C -21.6 ± 0.2‰). In seagrass-vegetated sediments, signatures ranged between δ¹⁵N 4.1 ± 0.2 (Grossenbrode) to 7.3 ± 0.3 (Falckenstein), and δ¹³C -23.1 ± 0.4 (Falckenstein) to -19.4 ± 0.5 (Wackerballig) ( Figure 4B ). δ¹³C and δ¹⁵N signatures were also homogenous across core depth; top sediments (0 to 5 cm) averaged across all seagrass cores were ¹⁵N depleted by 0.48‰ and ¹³C enriched by 0.68‰ compared to the 10 to 15 cm and 15 to 20 cm sections.

Averaged across all sites, visible organic material such as invertebrates and seagrass root, shoot, rhizome material on average contributed to 12% of the average total C_org in seagrass-vegetated sediments and 9% in unvegetated sediments. Visible organics were responsible for 65% of the C_org in unvegetated sediments in Hasselfelde. Results from the five-source two-biotracer (δ¹³C, δ¹⁵N) mixed model of the remaining organic fraction (the non-visible fraction, called SOC) in the sediments of Falckenstein (SOC fraction = 94%), Wackerballig (SOC = 56%), Grossenbrode (SOC = 78%), Graswarder (SOC = 55%), showed that phytoplankton (24%), P. littoralis (18%) and other macroalgae (22%), made the largest contribution to C_org overall, not seagrass biomass or its epiphytes, and a combination of two of these three dominant carbon sources could be seen at each site ( Figure 4A ). Sampling location accounted for most of the variation in the mixing model (50^th percentile σ = 6.311), while vegetation coverage (seagrass-vegetated vs unvegetated) had a negligible effect (50^th percentile σ = 0.105).

Radiocarbon dating of wood pieces found in Sierksdorf cores revealed no age-depth correlation within a single core, but all dates fell in two well separated time intervals, averaging at 5,806 years BP and 5,095 years BP (before present), i.e. a hiatus of approx. 700 years between these time intervals ( Table 3 ).

The percent fine sediment fraction of seagrass-vegetated sediments varied greatly between sites, from 0.09 ± 0.04% (Gahlkow) to 7 ± 3% (Goehren). All unvegetated sediments and 17 of 20 seagrass-vegetated sites had a low fine sediment fraction (<3%) ( Table 4 ). The unvegetated sediments of Orth and Seebar had the highest fine sediment fraction (2.5 and 3%), and of the seagrass-vegetated sites Glowe, Goehren, Niendorf, Orth, Sierksdorf had 2.7 to 7% fine sediments. Seagrass-vegetated sediments contained similar or more (by 1-31 times) fine-grained sediments than adjacent unvegetated sublocations, except in Graswarder, Gluecksburg, Seebar, where unvegetated sublocations contained 2-10 times more fine-grained sediments than their vegetated counterparts. For seagrass complexity, Maasholm exhibited the lowest (92 ± 21 m/m²) and Wackerballig the highest (539 ± 132 m/m²) average seagrass complexity of all sites (including sparse and dense seagrass sublocations) ( Table 4 ). Dense seagrass sublocations were one to five times more complex than sparse seagrass sublocations, with the greatest difference observed in Teichhof and Aschau, and smallest difference seen in Sierksdorf and Grossenbrode. Average MOV varied six-fold, ranging from 0.155 ± 0.007 m/s (Niendorf) and 0.91 ± 0.03 m/s (Teichhof) between sties, but some (Falshoeft lighthouse, Heidkate, Kellenhusen, Teichhof, Wackerballig) experienced strong peak MOVs (> 3 m/s) ( Table 4 ).

Six alternative candidate models were deemed statistically indistinguishable (ΔAICc < 2) from each other (summarized in Table 5A ). Meaning, there was no strong support for one particular model. In the averaged model, seawater depth, seagrass complexity, and fine sediment fraction were similarly important in predicting C_org within seagrass-vegetated sediments (RVI 1.00 for each, Table 5B ), and they had a significant negative (seawater depth) or positive (fine sediment fraction, seagrass complexity) effect on C_org stocks (see ‘Estimate’ in Table 5B ). Avg. MOV had a weaker (by approx. 2.5 times) and insignificant effect in predicting C_org. There were no significant interactive effects between variables, except for that of seawater depth and seagrass complexity, but it had a weak (by approx. 5 times) predictive effect on C_org stocks.

C_org stocks in the top 25 cm of seagrass-vegetated sediments in the German Baltic Sea were high (average 1,920 ± 402 g C/m²), and richer than adjacent unvegetated sediments, but differed widely between and even within sites. Regional heterogeneity observed here was comparable to previous regional evaluations of blue carbon in Z. marina meadows (e.g. Röhr et al., 2016; Prentice et al., 2019), but local heterogeneity was much greater, with previous measurements reporting <10 times more C_org content in seagrass-vegetated vs unvegetated sediments ( Table 6 ), or similar carbon content between these sublocations (Prentice et al., 2019; Prentice et al., 2020; Mazarrasa et al., 2021; Krause et al., 2022). Overall, average C_org stocks in the German Baltic Sea fit within the global range reported for those of Z. marina (318 ± 10 to 26,523 ± 667 g C/m²), but more closely resembled Z. marina C_org stocks in the Mediterranean Sea (8,793 + 2,248 g C/m²), Funen area of Denmark, and Skagerrak coast of Sweden, than C_org stocks in other parts of the Baltic Sea ( Table 6 ). In fact, German blue carbon stocks were on average richer (by three times) than the Baltic Sea average (578 ± 43 g C/m², excluding Germany, Röhr et al., 2018). A similar geographical trend was also reported for sites along the coast of Sweden, where percent C_org in top sediments were 10 to 25 times higher along the west (in Skagerrak) vs south and east coasts (in the Baltic Sea) of Sweden (Jephson et al., 2008). Consistent with the present study, wide regional variation in C_org stocks were observed in eastern Jutland of Denmark, which included “carbon hot spots” (stocks as high as 26,523 ± 667 g C/m² in Thurøbund) that were comparable to those observed in Germany (10,577 ± 6,178 g C/m²). Two biophysical parameters were thought to be responsible for the C_org hotspot found in Thurøbund: low wave exposure and high seagrass productivity (420 ± 98 shoots/m²) – seagrass densities and exposures that were similar in magnitude to those measured in the sheltered bay of Orth (467 ± 98 shoots/m²) where the third highest C_org was found in Germany. It is worthwhile contemplating explanations for the higher C_org stocks reported in Skagerrak-Kattegat and southwestern Baltic Sea relative to other parts of this basin ( Table 6 ). The rich geological history that shaped the Baltic Sea, including northward retreat of the Scandinavian ice sheet that caused land uplift in the south and changed coastal landscapes as a consequence of sea-level rise during the Littorina Transgression, but also shaped the geology that enhanced its ability to sequester C_org. For instance, hard-bottom and rocky shores dominate in the northern coasts, whereas till material and sandy/muddy beaches are common in the south (Schiewer, 2008).

Our results suggest that the C_org accumulating in seagrass meadows in the German Baltic Sea is primarily (88%) originating from allochthonous sources, material originating from outside the seagrass meadow, thus, most of the C_org making up stocks here are imported from outside the boundaries of the meadow. Of the five C_org sources tested, this material was predominantly derived from a combination of phytoplankton, drift algae P. littoralis, and other macroalgae. It must be noted that the determination of the source materials in the SOC fraction was derived from stable isotope analyses of materials from sites in Schleswig-Holstein only, and so the source material for Mecklenburg-Vorpommern may differ from these. Nonetheless, these findings are consistent with those of neighboring Baltic Sea nations, like Finland, where phytoplankton material was the primary source of C_org (43 – 86%), while seagrass made a relatively small contribution to the overall sediment C_org pool (1.5 – 32%) (Röhr et al., 2016), and along the northwest Pacific coast, Z. marina meadows primarily sequestered allochthonous C_org, originating from plankton, terrestrial, and kelp sources rather than seagrass material (Prentice et al., 2019; Krause et al., 2022). The latter contributed less than 25.3% to the C_org pool. However, in Denmark and Poland, seagrass biomass was found to be the biggest contributor to sediment C_org pools (13 – 81%, Röhr et al., 2016; 40 – 45%, Jankowska et al., 2016). A country-wide analysis of blue carbon stocks in Australia revealed that coastal meadows in temperate regions were dominated by autochthonous C_org sources (72% was derived from seagrass), whereas allochthonous material prevailed in tropical meadows (64%) (Mazarrasa et al., 2021). In a global analysis (across 88 locations and multiple seagrass species) of seagrass blue carbon stocks, seagrass biomass accounted for approximately 50% of the C_org stocks, while the remainder originated from allochthonous sources (Kennedy et al., 2010).

In the mixing model, sampling location (not vegetation coverage) was the main driver of the variation in C_org sources (for the non-visible fraction of C_org; SOC). Coastal landscape and differences in inputs to marine foodwebs, including currents, upwelling, as well as point (e.g. sewage outfall, rivers) and diffuse (e.g. atmospheric, rainwater runoff) sources of nutrients may help explain some of the dissimilarity observed between locations. For example, phytoplankton contribution was abundant in the seagrass meadow in Graswarder, but absent in Wackerballig. The former is situated near a point source contamination of sewage outflow (Schubert et al., 2013), and the latter is adjacent to a large marine national park. Phytoplankton composition and abundance are often used as bioindicators of enhanced nutrient concentrations (2000/60/EC, EU, 2000). For the drift algae P. littoralis, large mats can become entrapped in narrow fjords (Falckenstein) and bays (Wackerballig), whereas they may be carried away from sites along open coasts (Grasswarder, Grossenbrode), which fits with the pattern seen in the sediment sources of the present study.

While C_org content was consistently higher in seagrass-vegetated vs unvegetated sediments, two cases had similar (Gluecksburg and Hasselfelde) or more (Gohren and Seebar) C_org in their unvegetated sediments. A similar phenomenon was previously reported for Z. marina meadows and other temperate seagrass species (see Mazarrasa et al., 2017a; Mazarrasa et al., 2017b; Prentice et al., 2019; Prentice et al., 2020; Mazarrasa et al., 2021; Krause et al., 2022). Some of this discrepancy was attributed to the export of organic material originating from seagrass habitats to adjacent unvegetated sediments (Duarte and Krause-Jensen, 2017). A similar spillover effect of C_org from adjacent seagrass meadows may be contributing to the high C_org accumulating in nearby unvegetated sediments in Hasselfelde, where visible C_org (not SOC, but also not seagrass biomass) were the main source of C_org in unvegetated sublocations. Also, there is a dense community (several cm thick) of bivalve Cerastoderma edule populated unvegetated sediments at this site (but not in seagrass-vegetated sediments) and their soft tissue would have contributed to C_org sinks (or CO₂ sinks) measured in these sediments. Autochthonous export is not a valid justification for Gohren, Seebar and Gluecksburg as bare sediments were dominated by non-visible C_org (100% SOC) and the latter two had a sludge-like consistency with the top (Seebar) and sixth (Gluecksburg) highest fine (muddy) particle fraction of all sites and sublocations in the study. Heavy anthropogenic pressure at these sites may be contributing to the excess nutrient inputs that would elevate baseline C_org in the sediments (Nixon, 1995; Short and Burdick, 1996; Bowen and Valiela, 2001; Nedwell et al., 2002). Also, Seebar was not populated by seagrass only 8 years prior to the study, which could have contributed to these findings as well (pers. obs. P. Schubert). It must be noted that visible seagrass roots and rhizomes were omitted from this analysis, which may lead to a lower contribution of autochthonous organic carbon sources, as observed in the present study.

Unexpected large amounts of well-preserved wood pieces were found in one location: Sierksdorf, Luebeck Bay. Radiocarbon dating of this wood suggests that it was deposited here during two distinct time intervals, averaging at 5,806 BP and 5,095 BP, that coincide well with the second phase of the Littorina Transgression (dated approx. 6,000 to 3,800 BP) following the last deglaciation (Schmölcke et al., 2006; Kostecki et al., 2021). In this period, rapid changes occurred especially in the southwestern part of the Baltic Sea, where fast sea level rise combined with a flat landscape led to the widespread die off of forests situated along the coast. These events caused the demise of coastal alder woodlands that subsequently produced peat containing abundant alder-wood remains (Schmölcke et al., 2006). Indeed, the lack of age (¹⁴C) correlation across core depth is typical of natural wood pieces (rather than man-made artefacts) because younger wood can be found in deeper parts of the core due to the roots of the living tree. The two distinct time intervals occurring 700 years apart may suggest a recolonization event took place at this site. Interestingly, the two-time intervals are situated on either side of a brief cooling period (the Subboreal period), which started 5,650 BP and caused widespread natural environmental changes in the southwestern part of the Baltic Sea (Schmölcke et al., 2006). Note that while only one site was radiocarbon dated, we report old wood co-occurring with seagrass meadows in other parts of Schleswig-Holestein, Germany (see Table 1 ).

The oldest known submarine peatlands globally are dated 5,616 ± 46 years BP and correspond to the thick matte formed by Posidonia oceanica, a long-lived Mediterranean seagrass, which constitute deep and significant C_org stocks (Lo Iacono et al., 2008). Z. marina does not form similarly thick and old sedimentary deposits, but our investigation confirms that it too can store millennia of carbon by protecting former (and now submerged) forested peatlands, acting as a protective covering for these rich carbon deposits (Krause-Jensen et al., 2019). Even well-preserved prehistoric settlements have been previously discovered beneath Z. marina meadows in Denmark (Fischer, 2011; Andersen, 2013; Pedersen et al., 2017) and Germany (Goldhammer and Hartz, 2017). Submerged wood artefacts are also present in the Baltic Sea coast of Germany, including a site in Neustadt near Sierksdorf (Klooß, 2014), and they too overlap with seagrass habitats. Recognition that some seagrass meadows in Germany are sitting atop ancient terrestrial carbon deposits that are potentially several meters thick has important consequences for avoiding the re-emission of CO₂ (via C_org) stored on a millennial timescale. However, our understanding of submerged coastlines is incomplete, so it is likely that many more submarine peatlands, like that found in Sierksdorf, await discovery near the German Baltic Sea coast.

A fourth objective of the study was to use the identified relationships between seagrass meadow attributes and environmental parameters on C_org to extrapolate stocks for the entire German Baltic Sea region. This regional variation in C_org was best explained by seawater depth, seagrass complexity, and the fraction of fine particles in the sediment, not average MOV or interactions between parameters. The latter two had a positive effect on the C_org content in the sediment, whereas stocks decreased with seawater depth, suggesting that the highest stocks were found in shallow locations with high seagrass complexity and the ability to accumulate fine-grained particles, regardless of MOV at the seafloor. However, these parameters do not demonstrate a clear ability to predict the regional distribution of C_org ( Figure 5 ). A combination of these three parameters, along with water motion, is consistently found to be the main biotic and abiotic driver of C_org in the literature, but its influence varies widely across regions and seems principally coupled to local hydrodynamic regimes. Because sediment erosion and detritus export rates are typically heavily influenced by water motion (e.g. lower wave height and exposure, fetch, and currents), sediment C_org content is typically lower in dynamic systems compared to more static ones (e.g. Samper-Villarreal et al., 2016; Mazarrasa et al., 2017a; Prentice et al., 2019; Novak et al., 2020; Mazarrasa et al., 2021). But in the absence of strong water movement, like in the German Baltic Sea where tides are negligible, water currents are weak, and wave heights are low, other aspects of the systems can shape C_org stocks. The lack of a trend observed in our data in spite of sampling extremes in hydrodynamics of our system, is a testament to this hypothesis, whereby the sites that experienced weakest (Orth, Maasholm, Sierksdorf) maximum current speeds at the seafloor both had the highest C_org pools of all sites examined in this study. However, sites with the strongest (Teichhof, Wackerballig) currents did not have low C_org pools, in fact they had the 10^th and 11^th highest C_org of all sites. Reduced hydrodynamics may also mean enhanced particle trapping from the water column, which may help explain the mainly allochthonous provenance of the C_org and overall high stocks in our study.

The positive relationship observed between the amount of fine particles in the sediment and C_org in seagrass-vegetated sublocations is not surprising as it is well known that more C_org is associated with finer mineral particles in soils and sediments (Calvert et al., 1995; Lin et al., 2002). These fine particles nurture anoxic conditions in the uppermost layer of the soil that protect organic particles from remineralization (Schrameyer et al., 2018; Brodersen et al., 2019). Our findings are consistent with previous findings for Z. marina meadows in other than the German parts of the Baltic Sea (e.g. Sweden, Dahl et al., 2016; Denmark and Finland, Röhr et al., 2016) and elsewhere (Z. marina outside Baltic Sea: Miyajima et al., 2015; Miyajima et al., 2017; Prentice et al., 2019; Krause et al., 2022; multiple species, including Z. marina: Kennedy et al., 2010; P. oceanica: Gacia et al., 2002; Hendriks et al., 2008). For example, in Denmark and Finland, more than 40% of the variation in C_org between sites could be explained by sediment characteristics, including the fine sediment content.

The thick canopy of seagrass leaves is known to effectively intercept particles in the water column, as well as decrease sediment erosion and seagrass detritus export (Ward et al., 1984; Fonseca and Cahalan, 1992; Madsen et al., 2001; Christianen et al., 2013). While seagrass complexity is a strong positive predictor of the regional differences in C_org in our study, and that of others (Jankowska et al., 2016; Samper-Villarreal et al., 2016; Serrano et al., 2016), this was not the case for Z. marina meadows off the coast of British Columbia, in Canada, where no clear relationship (negative or positive) was observed between the two parameters (Prentice et al., 2019). Here, water motion was the strongest predictor of C_org.

The dampening effect of seawater depth on surface water motion is correlated to the trend of increasing C_org content with deeper depths (Lavery et al., 2013; Mazarrasa et al., 2017a). While we also observed a dampened effect of water currents (but at the seafloor, rather than surface waters) by increasing seawater depth (ρ = -0.43), an opposite relationship was seen between C_org content and depth. A fourfold decrease in C_org was also observed with increasing seawater depth (from 2–4 m to 6–8 m) in P. oceanica seagrass, in the Mediterranean Sea (Serrano et al., 2014). Regarding the Baltic Sea, seagrass was historically present at deeper depths (e.g. observed beyond 8 m in the past ten years), but today it is rarely observed beyond 5 m seawater depth in Germany. The realization of the Helsinki Commission’s Baltic Sea Action Plan (BSAP) goals, which implicate a considerable nutrient abatement, would result in seagrass expansion into deeper waters, as was shown by Bobsien et al., 2021 for the German part of the Baltic Sea. This expansion would lead to an enhancement of the CO₂ storage potential by these habitats, but it is important to consider that this potential is diminished at deeper depths (relative shallow depths) – an important consideration for carbon accounting when including seagrass blue carbon contributions to the German national CO₂ budget.

Our measurements confirm that seagrass meadows in the Baltic Sea coast of Germany store a large C_org pool. The high spatial heterogeneity seen across the region warrant site-specific investigations to obtain accurate estimates of blue carbon. However, localities with high seagrass complexity, high fine sediment fraction, and low seawater depth could help select localities with more favorable C_org accumulation potential. An unexpected and significant relic terrestrial C_org pool was found beneath the seagrass meadow in Sierksdorf (Luebeck Bay), and also confirmed in other locations (see Table 1 ). It is likely that many more submarine peatlands await discovery along the southwestern Baltic Sea region, and that they hold millennial timescale C_org deposits similar to those found in Germany.

Based on a conservative scaling up of measurements (integrated to 25 cm sediment depth), collectively (total of approx. 285 km², Schubert et al., 2015) seagrass meadows in the German Baltic Sea are preventing 2.01 Mt of future CO₂ emissions from being released into the atmosphere. However, it must be noted that in some locations (Gelting Bay, Teichhof) the sediment has been eroded to the marl layer, which seagrass roots cannot penetrate, while in other locations the sediment is known to extend to 7 m sediment thickness below the seafloor, e.g. within the inner areas of Mecklenburg Bay (Lemke, 1998), such as Sierksdorf, Niendorf, Kellenhusen, Grossenbrode in the present study.

Because C_org is dependent on high seagrass complexity, accumulation of blue carbon in Germany may be contingent on healthy seagrass habitats. Furthermore, loss of these habitats will have negative consequences for the German remaining CO₂ budget because the C_org stored beneath meadows may be rereleased into the water column and later to the atmosphere. Their loss would also impact their many co-benefits (see Heckwolf et al., 2021). Given the pressing need to offset and prevent future CO₂ emissions, more stringent and concerted efforts are urgently needed to enhance the C_org storage potential (via habitat restoration and improving growing conditions) and prevent further degradation (via conservation) of seagrass habitats along the Baltic Sea coast of Germany. Our study provides urgently needed knowledge and constitutes a further incentive to enhance efforts in protecting existing seagrass meadows in Germany, and restore them where natural recolonization is likely slow, such as in enclosed embayments or areas that have seen a loss in seagrass habitat, especially where the distance to the next seagrass-vegetated site is high.