Landscape- and site-scale spatial variability of blue carbon stocks and fluxes in tropical seagrass meadows

Seagrass meadows are emerging natural climate solutions for climate change mitigation through their high potential for organic carbon sequestration and storage, also known as blue carbon. However, the variability in current blue carbon stock and flux estimates is high, particularly at landscape scale. This knowledge gap highlights the need for evaluating blue carbon at spatial scales that are both locally robust and globally relevant. We quantified the magnitude of variability in blue carbon stocks and fluxes in tropical intertidal seagrass meadows at the site and landscape scales. We sampled within and across six intertidal seagrass meadows representing three geomorphic settings, including reef-associated settings dominated by Cymodocea spp., estuaries Halophila spp. and lagoons dominated by Thalassia spp., across Singapore. Within and across these sites, we measured soil organic carbon (C_org) stocks and greenhouse gas fluxes using the static chamber method. Tropical intertidal seagrass meadows stored 27.5±19.9 Mg C_org ha⁻¹ (mean ± SD) in the top 0–25 cm of soil, which varied significantly within sites (min–max: 5.1–78.9 Mg C_org ha⁻¹, n=42). Plot-level mean fluxes averaged 27.5±25.3 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ of CO₂ (min–max: −18.6 to 82.4, n=56) and 0.5±18.2 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ of CH₄ (min–max: −66.4 to 57.5, n=56). Mixed-effects models indicated that geomorphic setting did not explain variation in 0–25 cm stocks or in CO₂ and CH₄ fluxes. Sensitivity analyses for 0–100 cm stocks showed a marginal geomorphology effect for measured-only totals (p=0.074) that was not present after extrapolation to 100 cm (p=0.816). In multivariate driver models, CO₂ flux was negatively associated with water temperature (p=0.011) and showed an overall geomorphology effect (p=0.027), while CH₄ drivers were not resolved (all p≥0.158. Plot-level modelling further showed that CO₂ decreased with increasing 0–25 cm C_org stock (p=0.044), whereas CH₄ showed no association with C_org stocks (p=0.469). Overall, these results highlight strong spatial heterogeneity in tropical intertidal seagrass carbon dynamics and underscore the need for spatially and temporally explicit sampling and modelling frameworks that reflect nested spatial structure.

1 Introduction

Seagrass meadows have a high capacity for storing organic carbon in their sediments and are widely recognised as blue carbon ecosystems, where conservation and restoration are used to protect and enhance these sediment carbon stocks as a nature-based solutions for climate-change mitigation (Unsworth et al., 2022; Duarte et al., 2025; do Amaral Camara Lima et al., 2023). Yet seagrasses globally face a high risk of degradation, which could lead to ∼1154 Tg CO₂ emissions with a social cost of USD 213 billion (Krause et al., 2025a). Therefore, there is a global interest in incorporating the conservation and restoration of seagrass meadows into climate change mitigation schemes and policies, such as a country's Nationally Determined Contribution (NDC) to the Paris Agreement and associated Nationally Appropriate Mitigation Actions (NAMAs) (Macreadie et al., 2021). However, policies and management decisions should be guided by a robust evidence base, including credible estimates of the extent, drivers and patterns of organic carbon storage (Dahl et al., 2025).

Although knowledge and understanding of the magnitude and drivers of variation in organic carbon (C_org) storage in seagrass meadows has greatly increased recently at the global scale (Gomis et al., 2025; Krause et al., 2025a), there are still challenges in quantifying blue carbon at site-specific or landscape scales, which is a mosaic of habitats with unique surrounding biophysical environment, hindering global efforts to establish them as reliable natural climate solutions (Williamson et al., 2022; Stankovic et al., 2023). Variability in seagrass C_org and carbon budgets is driven by biotic and abiotic factors at various spatial scales (Krause et al., 2025a; Gomis et al., 2025; Mazarrasa et al., 2021; Mazarrasa et al., 2023). At the intra-site scale, extrapolating standard 100 cm C_org stocks from shallower cores is a major caveat, potentially leading to 1.5- to 10-fold overestimation (Dahl et al., 2025; Stankovic et al., 2023; Krause et al., 2025a). Local or site-scale variability of C_org is generally attributed to grain size (Miyajima et al., 2015), species composition (Gomis et al., 2025), the contribution of sources of organic matter, i.e. autochthonous and allochthonous (Ricart et al., 2020; Serrano et al., 2019), water depth, and hydrodynamic energy, (Lavery et al., 2013; Ndhlovu et al., 2024; York et al., 2018). However, our understanding of how these factors scale to broader landscape contexts remains incomplete.

At larger landscape scales, the coastal geomorphic setting has emerged as a key control on seagrass C_org variability (Kennedy et al., 2022; Mazarrasa et al., 2021); however, it is not known how the geomorphic setting would influence GHG emissions from seagrasses. Estuarine meadows generally inhabit the intertidal interface and shallower depths, receive large inputs of terrigenous organic material, and experience low hydrodynamic energy, creating a depositional environment (Ricart et al., 2020), potentially facilitating high accumulation of allochthonous carbon (Capece et al., 2025). Sheltered coastal settings – such as lagoons or back-reef meadows – can maintain higher surface-soil C_org due to reduced sediment resuspension and wave energy (Guerra-Vargas et al., 2020; Alemu et al., 2022). However, the same anoxic conditions and high organic matter inputs in seagrasses that promote high C_org burial and sequestration can also favour the production and emission of greenhouse gases such as methane (CH₄) (Rosentreter et al., 2021). Seagrass meadows generally emit negligible to moderate CH₄ fluxes (Asplund et al., 2022; Burkholz et al., 2020), although emissions of the higher magnitudes sometimes reported are more commonly associated with patchy seagrass cover associated with disturbed settings than with mature natural meadows (Halim et al., 2024). Because CH₄ flux magnitude is linked with geomorphic setting and its surrounding environmental controls such as organic matter loading, fine-sediment (mud) content, and electron acceptor availability – estuarine and sheltered lagoonal meadows, where C_org stocks are potentially higher, may also sustain conditions more conducive to methanogenesis than exposed or reef-associated settings. We note, however, that robustly assessing the climatic significance of these fluxes requires sustained accounting over the ecosystem lifetime e.g., 20- to 100-year Global Warming Potential horizons rather than comparisons of instantaneous flux magnitudes (Neubauer and Megonigal, 2015). Kirwan et al. (2023) proposed that GHG emissions in coastal wetlands would be governed by local-scale factors such as plant productivity and electron acceptor availability and linked environmental controls compared to landscape-scale geomorphic processes, yet this hypothesis has not been explicitly tested across geomorphic settings in seagrass meadows.

Here, we explore the spatial variability in carbon stocks and fluxes across seagrass meadows of different geomorphic settings in Singapore. We hypothesised that geomorphic settings such as reef-associated settings, estuaries and lagoonal settings would shape soil organic carbon and GHG emissions in tropical seagrass meadows, and a higher soil C_org would be associated with higher GHG fluxes under these geomorphic types. Here, we use “landscape scale” to refer to a mosaic of seagrass meadows distributed across different geomorphic settings (reef-associated, lagoonal, and estuarine) within Singapore, where larger-scale drivers of stocks and fluxes such as temperature and precipitation are held constant. We ask (1) what is the magnitude of soil organic carbon stocks and GHG emissions (CO₂ and CH₄) of tropical seagrass meadows, (2) How do organic carbon and GHG fluxes vary across spatial scales encompassing local (within‐site, across sites) and landscape (geomorphic settings) levels, and (3) are patterns in organic carbon storage associated with GHG fluxes, and which environmental parameters define these associations?

2 Methodology

2.1 Site description

The study was conducted on six mature natural intertidal seagrass meadows (except Eagle Bay, a reclaimed lagoonal site) across the tropical island state of Singapore (Fig. 1). The climate of the area is humid tropical; therefore, temperature does not vary much throughout the year, with the highest average monthly temperature of 28.8 °C during May and the lowest average monthly temperature of 26.8 °C during December and January (Meteorological Service Singapore, 2024). Mean annual rainfall is 2113 mm (1991–2024), with two monsoon periods during June–September and December–early March separated by inter-monsoon seasons in late March–May and October–November.

The seagrass meadows in this study are found in the lower intertidal zone and experience a diurnal tidal cycle, such that they are only briefly exposed during low tide (<0.7 m, tidal range >2.5 m) and remain fully submerged during high tide. We selected six intertidal meadows representing three geomorphic classes (McKenzie et al., 2016): reef-associated (Cyrene Reef, Labrador Nature Reserve; dominated by Cymodocea spp.), lagoonal (Tanah Merah, Eagle Bay; dominated by Thalassia spp.), and estuarine (Chek Jawa, Changi Beach Park, dominated by Halophila spp., Fig. 1). Table S1 in the Supplement summarises seagrass species composition and cover by site.

Figure 1 Locations of study sites across estuarine, lagoonal, and reef-associated geomorphic settings in Singapore (orange circles: estuarine; blue diamonds: lagoonal; green triangles: reef-associated). Seagrass extent across Singapore is shown on the main map but is barely visible owing to its limited areal cover; inset maps therefore show seagrass meadow boundaries (Tan et al., 2022) at higher resolution for the two largest studied meadows, Chek Jawa and Cyrene Reef.

2.2 Field sampling design

We collected field data during low tide from October 2022 to May 2024 and targeted the same time of the day and season for consistency and to avoid seasonal variations affecting GHG fluxes. We selected the sampling months when Singapore's climate is generally consistent with an average monthly temperature of 27.5–28.8 °C and average monthly rainfall of 150–160 mm (Meteorological Service Singapore, 2024). At each site, 3–5 transects (30 m length treated as within-site sampling units) were laid parallel to shore at approximately similar tidal exposure and elevations, >5 m apart, positioned to encompass the meadow extent following visual assessment. Three 1 m × 1 m plots were established along each transect at 5, 15, and 25 m. Where transects could not capture site heterogeneity, additional random plots were sampled within the meadow to improve coverage.

2.3 Carbon stocks

To quantify sediment carbon stocks, we collected 9–15 replicate sediment cores per site using a gouge auger (5.08 cm internal diameter) to 1 m depth or refusal. Because the upper 20–50 cm typically exhibits the strongest gradients in carbon content (Howard et al., 2014) due to the influence of the vegetation rooting zone, the top 50 cm was sectioned at 5 cm increments. Each core was divided into 11 sections: 0–5, 5–10, 10–15, 15–20, 20–25, 25–30, 30–35, 35–40, 40–45, 45–50, and 50–100 cm. We used 5 cm increments in the upper 50 cm to capture steep near-surface gradients commonly observed in seagrass sediments (Howard et al., 2014). Below 50 cm, we used a coarser interval (50–100 cm) due to expected lower variation in carbon stocks at depth (Phang et al., 2015). Core compaction during extraction was minimal (<10 %). Sediment was transported to the laboratory and dried at 60 °C to a constant mass. Dry bulk density (g dry weight cm⁻³) was calculated as dry mass divided by the in-situ sample volume (Howard et al., 2014). A subsample was ground and homogenised, then fumigated overnight with 0.1 M HCl to remove inorganic carbon, dried, weighed, encapsulated in tin, and analysed for total organic carbon (C_org, %) using an organic carbon analyser. Carbon stocks (Mg C_org ha⁻¹) were calculated following Howard et al. (2014), integrating section-specific C_org, dry bulk density, and layer thickness. We recorded core refusal depth (cm) during sampling and treated refusal depth as an informative sediment property reflecting sediment deposit thickness. Of the 62 cores collected, 25 (40.3 %) reached 100 cm depth based on refusal depth. To ensure comparable depth coverage across sites, primary inference for 0–25 cm stocks were restricted to cores that reached ≥25 cm depth (“complete-to-25” cores; n=42). We report 0–100 cm stocks as secondary sensitivity estimates (measured-only and extrapolated-to-100 cm), clearly labelled. For cores shorter than 100 cm, we extrapolated linearly from the deepest measured interval to 100 cm. For sediment texture, dried subsamples were disaggregated, sieved at 2 mm to remove coarse fragments (shells, pebbles), dispersed in deionised water, and analysed with a laser diffraction particle size analyser (Horiba LA, 350). Sediment texture was characterised using the mud fraction (% mud) measured for each depth interval; for modelling we summarised mud fraction as the mean % mud within 0–25 cm at the core level (stock models) and at the plot level (flux driver models).

2.4 Sediment-air greenhouse gas fluxes

We used opaque (dark) static chambers to quantify sediment–air CO₂ and CH₄ fluxes during low tide exposure representing sediment respiration/production, rather than net ecosystem exchange that would include photosynthetic uptake. Chamber collars (PVC, internal diameter 15 cm) were inserted ∼5 cm into the sediment. We minimised disturbance by careful insertion and waiting for a 20 min equilibration period before T₀ sampling based on our observations from a previous study (Iram et al., 2022). Chambers were rigid PVC with an enclosed headspace height of 30 cm. Lids were sealed using rubber bands and inspected sealing surface carefully to rule out leaks. Measurements were conducted once for each site between 10:00 and 17:00 SST (Singapore Standard Time) to minimise diel variability. Within each plot, three chambers were deployed at random positions and sealed with rubber bands Headspace gas was sampled by syringe at 0, 20, 40, and 60 min and transferred to 15 mL pre-evacuated glass vials (Labco, High Wycombe, UK). Concentrations were measured in the laboratory using a cavity ring-down spectrometer (Picarro G220-i). Fluxes (surface-area-normalised) were calculated from the linear rate of change in headspace concentration over time, corrected for chamber volume, surface area, ambient temperature and pressure via the ideal gas law. Flux estimates were retained when linear fits met r²≥0.8 ($\mathrm{162}/\mathrm{180}$) and non-linear estimates were excluded.

2.5 Environmental parameters

At each plot, we measured surface sediment temperature (°C) at the top 5 cm depth using an Electronic 4-in-1 soil meter. Sediment redox potential (mV), pH, salinity (psu), surface water temperature (°C), conductivity (mS cm⁻¹) and dissolved oxygen (%) were measured using a multiparameter Hanna^® pH/ORP sensor (HI769819X) immediately adjacent to chambers during low tide, allowing readings to stabilise for 5 min. Air temperature (°C) inside each chamber was measured using a Hobo^® Pendant MX Tem. All sensors were calibrated following the manufacturer's guidelines before each field campaign.

2.6 Data Analysis

All analyses were conducted in R (R Core Team, 2026) using lme4 (Bates et al., 2015) and lmerTest (Kuznetsova et al., 2017) for linear mixed-effects models (LMMs) and emmeans (Lenth and Piaskowski, 2026) for estimated marginal means. The primary stock response was 0–25 cm plot-level measured-only C_org stock restricted to 0–25 cm cores. Geomorphic setting effects on 0–25 cm stocks were tested using: Stock${}_{\mathrm{0}-\mathrm{25}}\sim$ geomorphology + (1 | site/transect).

Fluxes were analysed on the raw scale using Gaussian LMMs; negative fluxes were retained. CO₂ geomorphology-only models used plots nested within transects nested within sites (1 | site/transect/plot). For CH₄, the fully nested model was singular; therefore, the geomorphology-only inference used a stable structure (1 | site/plot). Fixed effects were tested using Type III tests with Satterthwaite denominator degrees of freedom (lmerTest).

To assess environmental controls on fluxes, we fit multivariate driver models including standardised (z-scored) predictors of redox, salinity, pH, water temperature, dissolved oxygen, and mud fraction (mean mud % in 0–25 cm), plus geomorphology. Conductivity was treated as redundant with salinity and was not included simultaneously. Pairwise comparisons among geomorphic classes and among sites were computed using Tukey-adjusted contrasts from estimated marginal means (emmeans); outputs are provided in the Supplement. We quantified how variability was distributed across spatial scales by extracting variance components from the fitted mixed models and expressing each component as a percentage of total variance (Table S2). We additionally explored an expanded within-site stock driver model including depth-resolved sediment properties (0–25 cm mud fraction and dry bulk density) and plot-level mean environmental covariates (salinity, redox, pH, temperature, dissolved oxygen), with predictors centred within site to separate within-site from between-site effects (Table S2). Table S2 provides full site-level variability outputs, including contrasts/EMMs with letters and N, refusal depth summaries, and variance partitioning.

2.7 Methodological limitations

We acknowledge methodological limitations due to logistical constraints. First, while we targeted high-resolution spatial replication within sites of ≤15 soil cores and 27 chambers per site, gas flux measurements are limited in temporal replication e.g., diel/seasonal variability and high-tide conditions and tidal range variability among sites were not sampled therefore the effect of driver responses which often operate on temporal scales such as salinity may constraint results. As such, we limit data reporting to per hr, not per day or per year. Second, below 50 cm depth we used a coarser core interval (50–100 cm) to prioritise replication, which may smooth fine-scale variability at depth and therefore interpretation of deep-layer estimates should be treated accordingly. However, based on prior measurements of seagrass meadows and other intertidal ecosystems in Singapore, we do not anticipate substantial variation in carbons stocks at lower depths (Phang et al., 2015). Third, chamber collars were inserted shortly before measurements due to access constraints during short low-tide windows; despite a 20 min equilibration period, some disturbance effects may remain. Fourth, we acknowledge that low-tide dark-chamber fluxes may differ under high-tide inundation due to altered oxidation/transport conditions; however, in this study, we are interested in sediment-air fluxes, rather than sediment-water fluxes.

3 Results

3.1 Soil carbon stocks and fluxes in tropical seagrass meadows

We standardised soil organic carbon (C_org) stocks to 0–25 cm as the primary reporting depth because 25 cm was the maximum common measured depth across sites. Across all cores, 0–25 cm C_org stocks averaged 27.5±19.9 Mg C_org ha⁻¹ (mean ± SD; range 5.1–78.9; n=42 cores; Fig. 2). Seagrasses had plot-level mean CO₂ fluxes of 27.5±25.7 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (Range: −18.6 to 82.4 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, n=54 plots, Fig. 3) and CH₄ fluxes of 0.5±18.2 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (Range: −66.4 to 57.5 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, n=54 plots; Fig. 4) for all measured sites.

Figure 2 Measured only 0–25 cm soil organic carbon stocks (Mg C_org ha⁻¹) across sites, restricted to the cores reaching ≥25 cm depth. Boxplots show medians and interquartile ranges; points show individual cores, coloured by geomorphology. Letters above sites indicate Tukey-adjusted differences in site estimated marginal means (emmeans; α=0.05); sites sharing a letter are not significantly different. N indicates the number of 0–25 cm cores per site.

Figure 3 Sediment–air CO₂ flux (mg m⁻² h⁻¹) across sites during low-tide exposure (dark chambers). Boxplots show chamber-level flux distributions; points show individual chambers, coloured by geomorphology. The y axis limits are symmetric (± max(|CO₂|) across all chambers) to aid cross-site comparison. Letters above sites indicate Tukey-adjusted differences in site estimated marginal means from the site mixed model (emmeans; α=0.05); sites sharing a letter are not significantly different. N indicates the number of plots contributing flux measurements at each site.

Figure 4 Sediment–air CH₄ flux ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) across sites during low-tide exposure (dark chambers). Boxplots show chamber-level flux distributions; points show individual chambers, coloured by geomorphology. The y axis limits are symmetric (± max(|CH₄|) across all chambers) to aid cross-site comparison. Letters above sites indicate Tukey-adjusted differences in site estimated marginal means from the site mixed model (emmeans; α=0.05); sites sharing a letter are not significantly different. N indicates the number of plots contributing flux measurements at each site.

3.2 Spatial variability in seagrass soil carbon stocks and fluxes

Mixed-effects models accounting for the nested sampling design showed no evidence of differences among geomorphic settings in 0–25 cm stocks (Fig. 2; Type III Satterthwaite test: $F(\mathrm{2},\mathrm{3.07})=\mathrm{0.246}$, p=0.796). CO₂ and CH₄ fluxes also did not differ among geomorphic settings in geomorphology-only models (CO₂, $F(\mathrm{2},\mathrm{3.00})=\mathrm{0.180}$, p=0.844; CH₄, $F(\mathrm{2},\mathrm{3.00})=\mathrm{0.844}$, p=0.513).

Sensitivity analyses for 0–100 cm stocks indicated a marginal geomorphology effect for measured-only totals ($F(\mathrm{2},\mathrm{2.65})=\mathrm{8.114}$, p=0.074), but no effect after extrapolating all cores to 100 cm ($F(\mathrm{2},\mathrm{2.96})=\mathrm{0.22}$, p=0.816). Refusal depth differed strongly among geomorphic settings (Table S2). Estuarine cores reached much greater depths (mean 97±15.3 cm, n=26) than lagoonal cores (mean 31.6±10.7 cm; n=18) and reef-associated cores (mean 18.7±8.2 cm; n=18). In a mixed-effects model with site as a random intercept, refusal depth differed significantly among geomorphic settings, with lagoonal (estimate −65 cm; p=0.003) and reef-associated (estimate $=-\mathrm{77.9}$ cm; p=0.002) cores being shallower than estuarine cores (Table S2). Overall, 25 of 62 cores (40.3 %) reached 100 cm depth (Table S2).

Variance partitioning indicated strong scale dependence in variability (Table S2): 0–25 cm stocks were dominated by site-level variance (72.30 %), with additional variance among transects within sites (7.61 %) and residual variance (20.09 %). CO₂variance was distributed across plot (16.96 %), transect (29.03 %), site (22.50 %) and residual (31.50 %) components, whereas CH₄ variance was concentrated at the plot scale (60.76 %) with smaller site (10.99 %) and residual (28.25 %) components.

3.3 Relationship between carbon stocks and fluxes and environmental parameters

Among tested sediment properties for 0–25 cm, mud % ranged from 2.5 % to 97.9 % and dry bulk density was 0.8 to 1.7 g cm⁻³ (see site and geomorphology level details in Table S3 and Fig. S1 in the Supplement). Across all sites, redox ranged from −91.6 to 156.7 mV, pH from 6.9 to 7.8, salinity from 22.4 to 31.3 psu, surface water temperature from 27.4 to 30.1 °C, conductivity from 33.4 to 48 mS cm⁻¹ and dissolved oxygen from 34.9 % to 97.7 % (see site and geomorphology level details in Table S4 and Fig. S2).

To explore drivers of site-dominated variability in 0–25 cm stocks, we fit an expanded plot-level mixed model including depth-resolved sediment properties (0–25 cm mud fraction and dry bulk density) and plot-mean environmental covariates (salinity, redox, pH, water temperature, dissolved oxygen). Predictors were centred within site (within-site effects) and transects were modelled as random intercepts nested within sites. None of the tested predictors were significant at α=0.05 (all p≥0.073); dissolved oxygen p=0.073, and salinity (p=0.099 were marginal (Table S2; StockDrv_Expanded sheets).

In the multivariate CO₂ driver model, geomorphology showed an overall effect (Type III Satterthwaite test: $F(\mathrm{2},\mathrm{16.42})=\mathrm{4.54}$, p=0.027). Adjusted EMMs showed that estuarine settings had significantly higher baseline sediment-surface CO₂ emissions (38.2±6.4 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) than lagoonal (22.1±5.1 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) and reef-associated settings (20.4±5.8 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; Tukey-adjusted p=0.038). These elevated estuarine emissions were associated with higher near-surface mud fractions and lower dissolved oxygen during low-tide exposure, while CO₂ fluxes across all settings declined significantly with increasing water temperature ($\mathit{\beta }=-\mathrm{12.4}$, p=0.011). Water temperature was a significant negative predictor (Type III Satterthwaite test: $F\mathrm{1},\mathrm{16.47})=\mathrm{8.12}$, p=0.011 while dissolved oxygen was marginal (p=0.081). Redox, salinity, pH, and mud fraction were not significant (all p>0.158). Tukey-adjusted pairwise contrasts are reported in Table S2. A sensitivity model removing the site random intercept produced an identical likelihood (likelihood-ratio test p=1.0 and slightly improved AIC (1467.866 vs. 1469.866), indicating robust fixed-effect inference (Table S2).

In the multivariate CH₄ driver model, no predictors were significant, including C_org (p=0.298), redox (p=0.703), salinity (p=0.262 and mud % (p=0.158). To evaluate stock–flux relationships at a consistent spatial scale, we compared plot-level 0–25 cm C_org stocks (Mg C_org ha⁻¹) with plot-mean GHG fluxes. In a linear mixed model with site as a random intercept, plot-mean CO₂ flux decreased significantly with increasing plot-level C_org stock ($\mathit{\beta }=-\mathrm{8.61}$, p=0.044; Fig. 5a), whereas plot-mean CH₄ flux showed no detectable association ($\mathit{\beta }=-\mathrm{2.25}$, p=0.469; Fig. 5b). Geomorphology contrasts were not significant for both CO₂ and CH₄ (all p≥0.399).

Figure 5 Relationships between measured 0–25 cm C_org stocks (plot-level) and plot-mean greenhouse gas fluxes. Panel (a) shows CO₂ versus 0–25 cm C_org stock; panel (b) shows CH₄ versus 0–25 cm C_org stock. Points represent plots coloured by geomorphology; lines show linear fits with 95 % confidence bands. Statistical inference is based on mixed-effects models that included geomorphology as a fixed effect and a random intercept for site.

4 Discussion

This study shows that tropical intertidal seagrass meadows in Singapore function as an important but highly spatially variable component of the coastal carbon cycle. This study quantifies near-surface (0–25 cm) soil organic carbon stocks and low-tide sediment–air CO₂ and CH₄ fluxes across six intertidal seagrass meadows in Singapore. Using a measured only 25 cm definition for near-surface stocks, geomorphology did not explain variation in 0–25 cm C_org stocks in mixed-effects models. Similarly, geomorphology-only models detected no setting differences in greenhouse gas fluxes (CO₂ and CH₄). However, in the multivariate CO₂ driver model that accounted for covarying environmental conditions (including temperature, oxygen and sediment texture), geomorphology showed an overall effect by structuring environmental gradients. Estuarine meadows promoted higher respiration than lagoonal and reef-associated settings through the combined influence of finer sediment texture (% mud) and lower surface-water oxygenation. This suggests that geomorphic settings may influence CO₂ flux indirectly through linked environmental gradients even when differences are not detectable in simpler setting-only tests. Our results highlight the necessity of spatially explicit carbon accounting that distinguishes between the drivers of storage and emissions to value these ecosystems in climate mitigation frameworks.

4.1 Soil carbon stocks and fluxes of tropical seagrass meadows in Singapore

Our soil organic carbon stocks for the top 15 and 30 cm fall within the ranges reported for local and global seagrass meadows (Table 1). Although our measured 100 cm soil organic carbon stocks are almost half of those reported for global seagrass meadows, they are consistent with those reported for the Halophila genus globally (Krause et al., 2025a). These findings situate intertidal seagrass meadows of Singapore as low to moderate carbon stores. The fluxes of CO₂ measured in this study are within the range of 5 to 920 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, reported for seagrasses in South and Southeast Asia (Stankovic et al., 2023; Zheng et al., 2023). An important observation from this study was the presence of negative fluxes from dark chambers, which had not been reported previously for seagrasses. However, this has been observed in other coastal ecosystems and linked to the presence of microbiota that can respire in the absence of light through novel metabolic pathway. This microbiota included acetogens, chemosynthetic microbes able to sequester CO₂ in dark conditions linked to chemoautotrophic pathways, in which chemolithoautotrophs and acetogens fix inorganic carbon without light via the Calvin-Benson-Bassham (CBB) and Wood–Ljungdahl pathway (Rodriguez et al., 2025). The potential presence of similar chemosynthetic microbes in the seagrass ecosystem suggests a yet unconstrained mechanism of CO₂ sequestration. Plot-mean CH₄ emissions averaged 0.5±18.2 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (range −66.4 to 57.5 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; n=54 plots; Fig. 4), within the lower range of reported values of sediment-water and sediment-air CH₄ fluxes for tropical and temperate seagrasses in the region and globally (Asplund et al., 2022; Henriksson et al., 2024; Zheng et al., 2023). The relatively low CH₄ emissions in our study likely reflect active methane consumption within sediment profiles. CH₄ produced in deeper anoxic zones is efficiently oxidized before reaching the atmosphere – by aerobic methanotrophs in surface oxic layers and by anaerobic methanotrophic archaea (ANME) operating in syntrophy with sulfate-reducing bacteria in anoxic sub-surface sediments (Orphan et al., 2001). These bacterial communities may consume up to 90 % of the CH₄ produced in marine sediments (Reeburgh, 2007). Our findings show that tropical seagrasses are sinks or negligible sources of CH₄.

Table 1 Comparison of seagrass soil carbon stocks with (15, 30 and 100 cm depth measured only unless indicated) with selected local and global estimates.

4.2 Scale-dependent variability in seagrass soil carbon stocks and fluxes

Across the primary geomorphology-only models, geomorphic setting was not a significant predictor of 0–25 cm stocks or of CO₂ and CH₄ fluxes. In contrast, geomorphology showed an overall effect in the multivariate CO₂ driver model after accounting for covariates, indicating that geomorphology may emerge through covarying environmental conditions rather than as a standalone setting effect. This indicates that, for the measured 0–25 cm stock response and the chamber fluxes analysed here, broad geomorphic classifications explain little additional variance beyond the site and within-site structure captured by the random effects. This is contrary to previous findings where geomorphic settings appeared as a strong predictor of 0–15 cm soil C_org stocks (Alemu et al., 2022), suggesting that geomorphic gradients may not be distinct enough due to the limited number of sites for each category or the influence of geomorphology may be overwhelmed by other factors. Refusal depth provides an informative indicator of sediment deposit thickness and potential capacity for deeper carbon storage. The strong geomorphology signal in refusal depth supports recommendations to report a common measured depth for primary inference in seagrass sediments and to treat deeper extrapolated estimates as sensitivity analyses. In our dataset, refusal depth patterns suggest that apparent differences in deeper (0–100 cm) stocks can partly reflect deposit thickness and sampling refusal rather than differences in near-surface carbon density alone.

Our findings highlight strong spatial variability in soil C_org stocks at the site level, which is consistent with differences in site history, disturbance, and terrestrial (allochthonous) inputs. Egale Bay, a reclaimed site with intensive recreational use, exhibited the lowest C_org stocks, whereas Tanah Merah – a sandy lagoon behind coastal revetments that receives substantial terrestrial inputs via a major drain – showed the highest C_org stocks.

High intra-site variability was likely due to site-specific patchiness in seagrass cover, species composition, and microtopography and species diversity. Seagrass species composition, cover, and canopy structure can influence sediment trapping, organic matter retention, and oxygen transport to sediments, thereby affecting both C_org storage (Kennedy et al., 2022; Krause et al., 2025b) and GHG production/oxidation. Although species traits were not explicitly modelled here, differences in meadow structure within and across sites could contribute to within-site and setting heterogeneity (Krause et al., 2025a) and should be incorporated in future trait-based blue carbon assessments. In our study, for instance, although Halophila spp. was dominant at both estuarine locations, Chek Jawa had higher diversity and C_org stocks than Changi Beach Park. Consistent with their C_org storage patterns, lagoonal and reef-associated meadows supported larger, high-biomass species (Thalassia spp. and Cymodocea spp.), contrasting with the smaller, low-biomass flora of estuarine sites (Table S4). This heterogeneity underscores the importance of adequate replication in sampling design to avoid misrepresentation of local stocks (Dahl et al., 2025). Generally, reef-associated and lagoonal meadows stored substantially more carbon than estuarine sites in the upper 15 cm soil profile (Table 1), consistent with the limited hydrodynamic energy environments of reef-associated and lagoonal meadows creating depositional conditions which favour carbon storage and retention (Kennedy et al., 2022; York et al., 2018). In contrast, estuarine meadows showed lower storage, similar to patterns reported in other estuarine meadows (Carruthers et al., 2007; Alemu et al., 2022). Despite large amounts of terrestrial runoff, site-specific exposure to strong hydrodynamic energy and sediment dynamics, such as resuspension of fine sediment, may reduce net C_org stocks in estuarine sites (Mazarrasa et al., 2023). In estuarine habitats, a balance of oceanic exchange and freshwater inflows creates a gradient in temperature, salinity and light that drives high variability in seagrass communities, which can influence intra-site C_org stocks in estuarine habitats (Carruthers et al., 2007; Alemu et al., 2022).

GHG fluxes exhibited substantial fine-scale heterogeneity. CO₂ variability was distributed across plot, transect, site, and residual components, while CH₄ variability was concentrated primarily at the plot scale, consistent with strong micro-scale patchiness and/or measurement-scale variability within sites. Similar to soil C_org storage, CO₂ and CH₄ fluxes showed limited variation across geomorphic settings at the landscape scale of Singapore, which exhibits a small spatial extent of 229.6 ha (Tan et al., 2022) thus represents an intermediate landscape scale. GHG fluxes were relatively consistent, suggesting that at intermediate landscape scales, GHG dynamics may be regulated by biogeochemical processes that are not influenced by geomorphology, geomorphic gradients may not be distinct enough at this scale, or the influence of geomorphology may be overwhelmed by other factors. Kirwan et al. (2023) proposed that fluxes are largely controlled by small-scale and short-term ecological and biogeochemical dynamics. Ecological processes such as species composition, plant productivity, electron acceptor availability and oxygen transport between soil and atmosphere can govern GHG fluxes in blue carbon ecosystems (Rosentreter et al., 2021). Our findings highlight that surface 25 cm carbon stock variability cannot be explained by geomorphic setting alone, in an intermediate landscape context and fluxes can appear similar across sites because they track production and deposition rather than long-term preservation dynamics.

4.3 Patterns and environmental controls on soil carbon stocks and fluxes

Our expanded core-level mixed model found that none of the tested sediment properties (mud fraction, bulk density) or plot-level mean environmental covariates explained site-dominated variability in 0–25 cm carbon stock. Although mud fraction is suggested as a proxy for intertidal seagrass meadows carbon stocks (Serrano et al., 2016), mud % was not a significant predictor of near-surface stocks or fluxes in our models, suggesting site-scale controls and unmeasured factors, e.g., hydrodynamic regime, deposit thickness, and species-specific traits (Kennedy et al., 2022) may be more important than within-site texture gradients in this system, rather than the within-site gradients captured by the measured physicochemical variables. In the multivariate CO₂ driver model, water temperature was the only significant predictor (negative effect), while mud %, redox, salinity, and pH were not significant and dissolved oxygen was marginal. The negative temperature relationship likely reflects temperature acting as a proxy for co-varying conditions (e.g., exposure/inundation state, timing within the low-tide window, or other unmeasured gradients) rather than a simple direct respiration-temperature response. For CH₄, none of the standardised predictors (including C_org %, mud %, redox, and salinity) was significant in the final driver model. These findings are consistent with previous studies where none of the potential drivers of vegetated coastal ecosystems, such as salinity, organic matter or biomass etc were associated with seagrass fluxes (Al-Haj and Fulweiler, 2020). However, we note that the absence of significant CH₄–salinity or CH₄–conductivity relationships must be interpreted within our sampling constraints. First, Singapore's polyhaline-to-euhaline seascape (22.4–31.3 psu) lacks the low-salinity conditions under which sulfate limitation typically enhances methanogenesis; within this narrow range, salinity exerts minimal control on CH₄ dynamics. Second, driver responses operate on temporal scales distinct from snapshot spatial measurements: low-tide static chambers capture instantaneous sediment–air exchange sensitive to short-term temperature and microtopography, whereas sediment carbon accumulation and microbial community structure integrate biogeochemical processes over seasonal to decadal scales. These scale mismatches may explain why classical macro-scale predictor–response relationships were absent at our intermediate landscape scale. Nonetheless, together with the variance partitioning (dominant plot-scale variance), this supports the interpretation that CH₄ fluxes in these intertidal meadows are highly heterogeneous at fine spatial scales and may require additional predictors e.g., microtopography, sulfate availability, plant traits, or episodic hydrologic forcing) and/or high-resolution temporal replication to resolve robust drivers. Overall, these findings reflect the nuanced interplay between environmental parameters and biogeochemical processes at multiple spatial scales, highlighting the importance of habitat-specific factors in controlling GHG emissions and organic carbon dynamics in tropical seagrass meadows and their implications for blue carbon budgets.

5 Conclusions

This study advances understanding of tropical seagrass carbon dynamics by quantifying soil organic carbon stocks and greenhouse gas fluxes across multiple spatial scales in Singapore. We show that 0–25 cm soil C_org stocks are highly variable but do not differ among geomorphic setting in mixed effects models. Geomorphology-only models similarly show no setting differences in CO₂ or CH₄ fluxes. In contrast, multivariate driver modelling indicates that CO₂ flux is associated with environmental gradients, with temperature as a significant predictor and an overall geomorphology effect after accounting for covariates, whereas CH₄ flux drivers were not resolved with the variables included. Refusal depth differed strongly among geomorphic settings, highlighting that sediment deposit thickness varies systematically across environments and supporting the use of common-depth primary inference alongside clearly labelled deeper sensitivity estimates. These findings underscore the importance of sampling designs and analyses that reflect nested spatial structure and of jointly considering carbon storage and greenhouse gas fluxes in robust coastal carbon accounting.

Our findings contribute critical baseline data from Southeast Asia – a region of high seagrass biodiversity but limited empirical data. They also provide a methodological framework for other small island states with comparable seascapes. Importantly, our study demonstrates that while tropical seagrasses can store substantial carbon, this ecosystem function is context-dependent, shaped by site-specific factors. Accounting for such variability is essential for robust integration of seagrasses into national and global climate strategies.