We thank the reviewer for their careful reading of our manuscript and for the constructive, detailed comments. We appreciate the reviewer’s suggestions to improve clarity around spatial scales, sampling design, environmental measurements, and the interpretation of greenhouse gas fluxes. We have utlined our point-by-point responses and planned changes in the attached PDF and we will revise the manuscript accordingly.

Thank you for the decision on the revisions and the opportunity to submit a revised manuscript. We have addressed all reviewer comments, which have substantially strengthened the manuscript. We are grateful to Referees #1 for their careful reading and constructive feedback. In response, we revised the manuscript to improve clarity on spatial scales, sampling design, environmental measurements, and interpretation of greenhouse gas fluxes. Key updates include:

• Basing primary inference for near-surface stocks on common measured depth across sites which was 0–25 cm stocks restricted to cores reaching ≥25 cm depth
• Redoing the statistics by using linear mixed-effects models (LMMs) and emmeans
• Adding depth-resolved mud fraction (% mud) to driver models and expanding the within-site stock driver model (mud fraction, dry bulk density, and plot-level environmental covariates)
• Updating plot-level stock–flux models
• Adding variance partitioning across spatial levels; and (vi) providing Supplementary workbooks (Table S1: study site information; Table S2: geomorphology contrasts, Type III ANOVAs, site variability, EMM letters/N, refusal depth summaries and models, and variance partitioning ((Referee 2)
• Figures 2-5 were regenerated to show distributions, sample sizes, and model-based letter groupings (Referee 1 & 2)
• Added seagrass species composition data in methods and context in discussion

Referee #1

1. Comment: Title misleading; comparison is among geomorphic settings and sites (two spatial scales), not just landscape. Suggest changing the title.

Response: We have revised the title to reflect both site- and geomorphic-scale comparisons and specify intertidal seagrass meadows.

Manuscript changes: Title updated (Title; Page 1, Lines 1–2).

• Old: “Landscape-scale spatial variability of blue carbon stocks and fluxes in tropical seagrass meadows”
• New: Landscape- and site-scale spatial variability of blue carbon stocks and greenhouse gas fluxes in tropical intertidal seagrass meadows

1. Comment: Abstract: name the species here, or total number of species involved in the case of multi-specific meadows.

Response: We have added species names and details in Methods and Supplementary Table S1.

Manuscript changes: Abstract + Methods + Supplement Table S1.

Revised manuscript text:

• Abstract; Page 1, Lines 15–17. We sampled within and across six intertidal seagrass meadows representing three geomorphic settings, including reef-associated settings dominated by Cymodocea spp., estuaries dominated by Halophila spp. and Cymodocea spp. and lagoons dominated by Thalassia spp. Page1, Line 15-18.
• Methods; Page 3, Lines 84–86. 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). Page 3, Line 87-89.

1. Comment: Add “in seagrasses”; reverse word order.

Response: We have added “in seagrasses” and reversed the word order as suggested.

Manuscript changes:

Revised manuscript text:  Introduction (Page 2, Line 33-35). Yet globally, seagrasses face a high risk of habitat loss and degradation, which could lead to ~ 1,154 Tg CO2 emissions with a social cost of $213 billion USD (Krause et al., 2025).

1. Comment: “Landscape” definition: better described as a mosaic of habitats.

Response: We have clarified “landscape scale” as a mosaic of meadows distributed across geomorphic settings within Singapore.

Manuscript changes: Introduction (Page 3, Lines 70–73). 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.

1. Comment: Add “C_org” formatting/clarity.

Response: We have standardised notation and ensured consistent use of C_org throughout.

Manuscript changes: Throughout the manuscript.

Revised manuscript text example. Introduction; Page 2. Lines 44-47. Variability in seagrass C_org and carbon budgets is driven by biotic and abiotic factors at various spatial scales (Krause et al., 2025; 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; Krause et al., 2025; Stankovic et al., 2023).

1. Comment: Methodology: no explanation how environmental variables were measured.

Response: We have added a dedicated subsection detailing instrumentation and units; and updated Results ranges and Supplement tables accordingly.

Manuscript changes: Methods (2.5 Environmental parameters) + Results (3.3) + Supplement Table S4 + Fig. S2.

• Revised manuscript text. Methods. Page 6, Lines 141-146. 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 minutes. 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.
• Results: Page 10, Lines 227-232. Among tested sediment properties, mud % ranged from 2.5 to 97.9% and dry bulk density was 0.8 to 1.7 g cm-3 (see site and geomorphology level details in Table S3 and Fig S1). 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).

1. Comment: Fig 1: add site labels; symbol count looks like 9 sites; avoid hiding text with zoom.

Response: We have redesigned Figure 1 with correct site labels and clearer symbology.

Manuscript changes: Figure 1 + caption revised.

Revised manuscript text: Page 4, Line 92-94. Revised caption. Extent of seagrass across Singapore (Tan et al., 2022) and study sites locations of intertidal seagrass meadows across estuarine, lagoonal and reef-associated geomorphic settings. Two sites with the largest seagrass extent among the studied sites are zoomed in.

1. Comment: Field sampling design: different exposure to air? Same depths across sites?

Response: We have clarified comparable tidal exposure/elevation where feasible and that flux measurements were during low-tide exposure.

Manuscript changes: Methods 2.2 (Page 4, Line 100-103). 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.

1. Comment: Transects as replicates? clarify.

Response: We have clarified the hierarchy and treated transects as within-site sampling units, not independent site replicates.

Manuscript changes: Methods 2.2 (Page 4, Line 94). 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.

1. Comment: Core slicing: 5 cm slices not standard; 50–100 cm length could be variable—justify design.

Response: We have justified finer resolution in the upper profile and acknowledged limitations of the coarse deep interval.

• Manuscript changes: Methods (Carbon stocks; Page 4–5, Lines 106–111). 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).
• Methods (Section 2.7; Page 7, Lines 172–173. 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.

1. Comment: Extrapolated to 1 m: how short? how many?

Response: We have reported refusal depth (cm), set measured-only 0–25 cm stocks (0-25cm; n=42) as primary inference, and presented 0–100 cm (measured-only and extrapolated-to-100 cm) as sensitivity analyses. We have also reported that 25/62 (40.3%) cores reached 100 cm and added refusal depth summaries/models in Table S2.

Manuscript changes:

• Methods 2.3. Page 5, Lines 117-123. Our primary carbon stock response was measured-only soil organic carbon stock in the 0–25 cm interval. To ensure comparable depth coverage across sites, primary inference for 0–25 cm stocks was 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.
• Results 3.2 Page 9, Lines 214-221. Sensitivity analyses for 0–100 cm stocks indicated a marginal geomorphology effect for measured-only totals (F (2,2.65) = 8.114, p = 0.074), but no effect after extrapolating all cores to 100 cm (F (2,2.96) = 0.22, p = 0.816). Refusal depth differed strongly among geomorphic settings (Supplement 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 = -77.9 cm; p = 0.002) cores being shallower than estuarine cores (Supplement Table S2). Overall, 25 of 62 cores (40.3%) reached 100 cm depth, and 25 of 62 cores (40.3%) reached ≥80 cm depth.

1. Comment: Low-tide sampling: add how results would differ in high tide (literature or limitations).

Response: While it is common to measure soil-air vertical fluxes at low tide, we added a limitation noting that low-tide dark-chamber fluxes may differ under high-tide inundation due to altered oxidation/transport conditions.

Manuscript changes: Methods 2.7 (Page 7, Lines 177–189). We acknowledge that 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.

1. Comment: Only dark chambers? What about photosynthesis?

Response: We have clarified that dark chambers represent exposed sediment-community respiration/production rather than net ecosystem exchange.

Manuscript changes: Methods 2.4 (Page 5, Lines 129–130). We used opaque (dark) static chambers to quantify sediment–air CO2 and CH4 fluxes during low tide exposure, representing sediment respiration/production, rather than net ecosystem exchange that would include photosynthetic uptake.

1. Comment: What time of day were chamber measurements done?

Response: We have reported the measurement window (10:00–17:00) and rationale.

Manuscript changes: Methods 2.4 (Fluxes; Page 5, Line 1333-134).

Revised manuscript text: Measurements were conducted between 10:00 and 17:00 to minimise diel variability.

1. Comment: Collar insertion affects GHG; ideally install 1–2 weeks prior—was this done?

Response: We have stated collars were inserted shortly before measurement due to access/safety constraints; we used a 20-minute equilibration period based on our observation from previous studies (Iram et al., 2021) and list this as a limitation.

Manuscript changes: Methods 2.4 (Fluxes; Page 5, Lines 131–133).

Revised manuscript text: We minimised disturbance by careful insertion and waiting for a 20-minute equilibration period before T0 sampling based on our observations from previous study (Iram et al., 2021).

1. Comment: Chambers flexible or rigid? height? airtight checks? Please explain chamber design.

Response: We have expanded the chamber description (rigid PVC; headspace height 30 cm; sealing and leak checks).

Manuscript changes: Methods 2.4 (Fluxes; Page 5, Lines 133–135).

Revised manuscript text: Chambers were rigid PVC with an enclosed headspace height of 30 cm. Lids were sealed using rubber bands and checked for leaks by repeated ambient readings.

1. Comment: CH4 outliers: cleaned before analysis? plots raw vs transformed?

Response: We have retained negative fluxes and did not remove statistical outliers a priori; we applied objective QC (e.g., r2 threshold) and used Gaussian LMMs on the original scale. Figures show raw distributions.

Manuscript changes: Methods (Data analysis, Page 6, Line 152)

Revised manuscript text:  Methods (Data analysis). Fluxes were analysed on the raw scale using Gaussian LMMs; negative fluxes were retained.

1. Comment: Transformations: pick one; if multiple, specify; add a Supplement table.

Response: We did not use any transformation in our revised analysis. We have used Gaussian LMMs on the original scale (no transformation), with negative fluxes retained; QC is described.

Manuscript changes: Methods (Data analysis; transformation text removed, Page 6, Line 152).

Revised manuscript text: Methods (Data analysis). Fluxes were analysed on the raw scale using Gaussian LMMs; negative fluxes were retained.

1. Comment: Use GLM/GLMM rather than Kruskal–Wallis; same for correlations.

Response: We have rerun our statistics; primary inference uses Gaussian LMMs reflecting the nested hierarchy for stocks, fluxes, and drivers.

Manuscript changes: Methods (Data analysis Page 6, Lines 148–166)

Revised manuscript text: 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 measured-only C_org stock restricted to complete-to-25 cores. Geomorphic setting effects on 0–25 cm stocks were tested using: Stock0-25 ~ 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 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 standardized (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 (Supplement 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 environmental covariates (salinity, redox, pH, temperature, dissolved oxygen), with predictors centred within site to separate within-site from between-site effects (Supplement Table S2). Supplement Table S2 provides full site-level variability outputs, including contrasts/EMMs with letters and N, refusal depth summaries, and variance partitioning.  

1. Comment: Add references for R and packages.

Response: We have added citations for R and key packages.

Manuscript changes: Methods 2.6 (Data analysis; Page 6, Lines 148–149).

Revised manuscript text: 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.

1. Comment: Figures: add N; clarify plot types; add df/test type; remove/move redundant figures to Supplement.

Response: We have updated figures to nested boxplots with points, included N and model-based letters (Tukey-adjusted EMM comparisons), and provided additional outputs in Table S2.

Manuscript changes: Figs. 2–5 + captions updated; Pages 8, 9, 11; Supplement Table S2.

Revised manuscript text: Caption Figure 2. Measured-only 0–25 cm soil organic carbon stocks (Mg Corg ha−1) across sites, restricted to complete-to-25 cores (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 complete-to-25 cores per site.

1. Comment: Significant differences among sites/settings in C stocks? Why not analyse?

Response: We have analysed this in our revised methods; geomorphology was not significant for 0–25 cm stocks (0-25cm; F (2,3.07) = 0.246, p = 0.796).  Site-level differences are reported via EMM letters and N (Table S2).

Manuscript changes: Results 3.2 + Fig. 2 (Page 8, Lines 186-1190) + Table S2.

Revised manuscript text: Results 3.2. Page 9, Lines 210-211. 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 (2,3.07) = 0.246, p = 0.796).

1. Comment: Grain size: add % mud (<63 µm) or range; add results mean/range.

Response: We have included depth-resolved mud fraction (% mud), incorporated mud into driver models, and reported results accordingly.

Manuscript changes: Methods 2.3 and 2.6 + Results 3.3 + Table S3 + Fig. S1.

Revised manuscript text:

• Methods 2.3. Page 5, Lines 126-128. 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).
• Methods 2.7, Page 6, Lines 157-158. To assess environmental controls on fluxes, we fit multivariate driver models including standardized (z-scored) predictors of redox, salinity, pH, water temperature, dissolved oxygen, mud fraction (mean mud% in 0–25 cm) and plus geomorphology.
• Results 3.3. Page 9, Line 229-230. Among tested sediment properties, mud % ranged from 2.5 to 97.9% and dry bulk density was 0.8 to 1.7 g cm^-3 (see site and geomorphology level details in Table S3 and Fig S1).

1. Comment: Is any driver significant? If yes, state it.

Response: Yes, water temperature was significant CO₂ driver, and we have revised the text to state this clearly. We have reported significant predictors and model outputs for CO₂ and CO₄ driver models and stock–flux relationships.

Manuscript changes: Results 3.3 + Fig. 5 + Table S2.

Revised manuscript text (Page 10, Lines 232–248): To explore drivers of site-dominated variability in 0–25 cm stocks, we fit an expanded core-level mixed model including depth-resolved sediment properties (0–25 cm mud fraction and dry bulk density) and plot-level 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 (Supplement Table S2; StockDrv_Expanded sheets).

In the multivariate CO₂ driver model, geomorphology showed an overall effect (Type III Satterthwaite test: F (2,16.42) = 4.54, p = 0.027). Water temperature was a significant negative predictor (Type III Satterthwaite test: F 1,16.47) = 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 Supplement 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 (Supplement 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). At the plot level (plot-summed 0–25 cm stock vs plot-mean flux) with site as a random intercept, CO₂ flux decreased significantly with increasing plot-level stocks (β = -8.61, 95% CI = [-17.0, -0.26], p = 0.044: Fig. 5a). CH₄ flux showed no detectable association with 0–25 cm stock in the analogous plot-level model (β = -2.25, 95% CI = [8.57, 4.06], p = 0.469: Fig. 5b). Geomorphology contrasts were not significant for both CO2 and CH4 (all p ≥ 0.399).

1. Comment: Discussion short; add take-home message and move/add methodological limitations.

Response: We have revised the discussion to include a clearer take-home message that tropical intertidal seagrass meadows were highly spatially variable and geomorphology did not explain variation in 0–25 cm C_org stocks and greenhouse gas fluxes (CO₂ and CH₄) in mixed-effects models analysis. We have added a methodological limitations section in methods, covering aspects such as temporal replication, core slicing below 50 cm depth, chamber collar insertion and low-tide sediment-air fluxes.

Manuscript changes: Discussion + Limitations updated.

Revised manuscript text:

• (Page 12, Lines 255–261) This study shows that tropical intertidal seagrass meadows in Singapore function as 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 CO2 and CH4 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 Corg stocks in mixed-effects models. Similarly, geomorphology-only models detected no setting differences in greenhouse gas fluxes (CO2 and CH4). However, in the multivariate CO2 driver model that accounted for covarying environmental conditions (including temperature, oxygen and sediment texture), geomorphology showed an overall effect.
• Methodological limitations. (Page 7, Lines 169–189). We acknowledge methodological limitations due to logistical constraints. First, while we targeted high-resolution spatial replication within sites (≤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). 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-minute 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.

1. Comment: Add within-estuary variability/estuarine gradient discussion.

Response: We have added discussion of estuarine gradients and within-estuary variability as contributors to heterogeneity.

Manuscript changes: Discussion (Page 14, Lines 318–323).

Revised manuscript text: 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).

1. Comment: Was seagrass cover measured? species present? canopy differences?

Response: Yes, we have added our recorder data of species diversity and cover context from literature and framed species traits as plausible contributors to heterogeneity; trait data were not explicitly modelled.

Manuscript changes: Discussion (Page 14, Lines 305–313) + Supplement Table S1.

Revised manuscript text:

            Comment: 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 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., 2025) 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).

1. Intermediate scales unclear (m or km?)—specify.

Response: We have clarified this by adding a numeric landscape extent for Singapore seagrass area.

Manuscript changes: Discussion (Page 14, Lines 327–329).

Revised manuscript text: 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 seagrass spatial extent of 229.6 ha (Tan et al., 2022), thus representing an intermediate landscape scale.

1. Comment: Add discussion/limitations for high tide.

Response: Though it is common to measure water-air fluxes during high tide, we focused on the sediment-air flux interface in the present study and have explicitly mentioned this throughout the manuscript. Regardless, we have added a high-tide/inundation limitation and implications.

Manuscript changes: Methods 2.7 (Page 7, Lines 177–189).

Revised manuscript text: 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.

1. Comment: Statement speculative; rewrite.

Response: We have revised wording to avoid speculation and interpret associations cautiously.

Manuscript changes: Discussion (Page 14, Lines 327–332).

• Old text: Our findings highlight that carbon stock variability can be explained by geomorphic settings (estuarine vs barrier-lagoon systems, sediment supply, tidal range) even in an intermediate landscape context, while fluxes can appear similar across sites because they track production and deposition rather 290 than long-term preservation dynamics.
• Revised text: 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 seagrass spatial extent of 229.6 ha (Tan et al., 2022) thus, representing an intermediate landscape scale. GHG fluxes were relatively consistent, suggesting that at intermediate landscape scales, either 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.

1. Comment: Add references; clarify “sources of carbon.”

Response: We have added the reference and removed the vague statement in the revised discussion for CH₄ drivers.

Manuscript changes: Discussion 4.3 (Page 15, Lines 332–335).

Revised manuscript text: 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).

1. Comment: Seagrass biology missing; add ecological processes/literature.

Response: We have expanded the Discussion to integrate species traits/canopy structure and ecological context alongside biogeochemical drivers.

Manuscript changes: Discussion expanded 4.2 (Page 14, Lines 305–313).

Revised manuscript text: 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 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., 2025) 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).

Referee #2

Thank you for the decision on the revisions and the opportunity to submit a revised manuscript. We have completed all reviewer comments as requested, and they have improved our manuscript considerably. We have made some substantial edits, and a summary of the biggest changes we made include:

• Basing primary inference for near-surface stocks on common measured depth across sites which was 0–25 cm stocks restricted to cores reaching ≥25 cm depth
• Redoing the statistics by using linear mixed-effects models (LMMs) and emmeans
• Reporting refusal depth (cm) as an informative result and adding 0–100 cm sensitivity analyses (measured-only and extrapolated-to-100 cm)
• Updating plot-level stock–flux models
• Adding variance partitioning across spatial levels; and (vi) providing Supplementary workbooks (Table S1: study site information; Table S2: geomorphology contrasts, Type III ANOVAs, site variability, EMM letters/N, refusal depth summaries and models, and variance partitioning
• Figures 2-5 were regenerated to show distributions, sample sizes, and model-based letter groupings
• Added seagrass species composition data in methods and context in discussion

1. Comment: (General) The manuscript by Iram et al. samples and estimates organic carbon (Corg) stocks in intertidal seagrass meadows in Singapore across estuarine, reef, and lagoonal settings. The study also includes CO2 and CH4 flux measurements using static chambers. Information on Southeast Asian seagrass ecosystems remains relatively scarce, particularly regarding GHG emissions, so these data will likely be well received by the community. The main objective of the manuscript is to assess spatial variability in carbon stocks and GHG fluxes across multiple scales. In my view, however, the way spatial variability has been assessed represents the main concern of this review, as it affects both the interpretation and visualization of the results. To evaluate variability across scales, the sampling structure and statistical treatment should reflect the study’s spatial hierarchy (e.g., geomorf. setting, site, plot, and replicate measurements). At present, the number of samples at each level and the nested design are unclear.  The sampling hierarchy could be better described in the Methods. A diagram or table summarizing the nested sampling design would be helpful for readers. For example, the authors could outline the spatial levels considered (e.g., geomorf setting, site, plot, and measurement) and the replication at each level, indicating how many sites were sampled within each setting (estuarine, reef, lagoon), how many plots were established per site, and how many sediment cores or chamber measurements were collected per plot. This would give some structure at how spatial variability was assessed.

Response: We have documented the spatial hierarchy and reflected it in nested LMMs and added site-level variability outputs (pairwise contrasts, EMM letters, N) and geomorphology contrasts/ANOVAs in methods, results and Supplement workbooks (Table S2 added).

Manuscript changes. Revised manuscript text:

• Methods 2.6 (Page 6, Line 149-155). The primary stock response was 0–25 cm measured-only C_org stock restricted to 0-25 cm cores. Geomorphic setting effects on 0–25 cm stocks were tested using: Stock_0-25 ~ 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).
• Results 3.2 (Page 9; Lines 210-226). 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 (2,3.07) = 0.246, p = 0.796). CO₂and CH₄ fluxes also did not differ among geomorphic settings in geomorphology-only models (CO₂, F (2,3.00) = 0.180, p = 0.844; CH₄, F (2,3.00) = 0.844, p = 0.513). Sensitivity analyses for 0–100 cm stocks indicated a marginal geomorphology effect for measured-only totals (F (2,2.65) = 8.114, p = 0.074), but no effect after extrapolating all cores to 100 cm (F (2,2.96) = 0.22, p = 0.816). Refusal depth differed strongly among geomorphic settings (Supplement 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 = -77.9 cm; p = 0.002) cores being shallower than estuarine cores (Supplement Table S2). Overall, 25 of 62 cores (40.3%) reached 100 cm depth, and 25 of 62 cores (40.3%) reached ≥80 cm depth (Supplement Table S2). Variance partitioning indicated strong scale dependence in variability (Supplement 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.

1. Comment: In the Results, the current data visualization approach, where all replicates are averaged per plot and per site to build a single depth profile, may not be the most effective way to represent spatial variability. To better convey variability across plots, sites, and settings, I suggest using boxplots within depth segments at each site, and facets to separate panels by setting to reduce clutter.

Response: We have revised Figure 2 using boxplot to preserve within-site distributions and avoid over-aggregation in the suggested boxplot style with facets by settings.

Revised manuscript text:

Figure 2 revised. Page 8, Lines 186–194. Caption: Figure 2. Measured only 0–25 cm soil organic carbon stocks (Mg C_org ha⁻¹) across sites, restricted to 0-25 cores (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.

1. Comment: The Corg stock bar plots in Fig. 3 could be more effectively represented as nested boxplots. I suggest displaying individual plot-level measurements as points overlaid on the boxplots (to improve transparency), grouped by site and colored by setting. This keeps variability within and among sites and settings, while preserving the nested structure of the data. Additionally, significant differences for stocks could be indicated using letters directly on the plot. The same approach could be applied to flux measurements. A single nested plot per variable (stocks or fluxes) should be sufficient, making the separate boxplots for each geomorphic setting (current panel B) unnecessary.

Response: We have replaced bar plots with nested boxplots + points, added letters and N, and removed redundant panels.

Manuscript changes: Figs. 2–5 revised; captions updated (Page 8–9, Lines 187–207).

Revised manuscript text:

• Updated caption Figure 2. Measured only 0–25 cm soil organic carbon stocks (Mg C_org ha⁻¹) across sites, restricted to 0-25 cores (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.
• Updated caption Figure 3. Sediment–air CO2 flux (mg m−2 h−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(|CO2|) across all 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.
• Updated caption Figure 4. Sediment–air CH4 flux (µg m−2 h−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(|CH4|) 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.

1. Comment: The results, as currently written, read somewhat redundant, because the same text structure is repeated for each variable, starting with the overall mean of all measurements, then reporting a percentage describing intra-site heterogeneity, followed by highlighting the sites with the highest and lowest values, noting differences between geomorphic settings, and mentioning vertical distribution patterns where relevant. While this is descriptive, a layer of analytical depth could be added. For example, the authors could quantify spatial variability using metrics such as the coefficient of variation and variance partitioning (e.g., nested ANOVA or mixed-effects models) to determine how much variability is explained at each spatial scale. Currently, the mixed-effects models are only used to test environmental controls on GHG fluxes and stocks, but they could also be leveraged to assess the nested spatial structure of the data.

Response: We have reduced redundancy by removing relevant text and figures. We have added variance partitioning from mixed models to add analytical depth as suggested by reviewer (reported in Table S2). We have also used LLM for assessing spatial structure of the data as described above in comment 1 response.

Manuscript changes: Results 3.2 updated; Table S2 variance partitioning sheets added.

Revised manuscript text: Page 10, Lines 222–226. Variance partitioning indicated strong scale dependence in variability (Supplement 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.

1. Comment: Spatial variability: unclear why measurements are lumped across plots/sites/settings.

Response: We have revised the results and figures to report distributions by site and avoid pooled summaries as primary inference.

Manuscript changes: Results 3.2 revised (Page 9–10, Lines 201–219); Table S2 provides contrasts/EMMs.

Revised manuscript text: 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 (2,3.07) = 0.246, p = 0.796). CO₂ and CH₄ fluxes also did not differ among geomorphic settings in geomorphology-only models (CO₂, F (2,3.00) = 0.180, p = 0.844; CH₄, F (2,3.00) = 0.844, p = 0.513).

Sensitivity analyses for 0–100 cm stocks indicated a marginal geomorphology effect for measured-only totals (F (2,2.65) = 8.114, p = 0.074), but no effect after extrapolating all cores to 100 cm (F (2,2.96) = 0.22, p = 0.816). Refusal depth differed strongly among geomorphic settings (Supplement 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 = -77.9 cm; p = 0.002) cores being shallower than estuarine cores (Supplement Table S2). Overall, 25 of 62 cores (40.3%) reached 100 cm depth, and 25 of 62 cores (40.3%) reached ≥80 cm depth (Supplement Table S2).

1. For each soil parameter (e.g., stocks, Corg content, Corg density, bulk density, grain size), the results are reported as mean ± SD or SE (it is unclear which) and the range of all measurements, lumped together across plots, sites, and settings, at the start of each Results section for the parameter of interest. Several points here:

• Depth basis: Up to which depth were these means calculated? Were they based on a common depth across all cores? Do the reported values include extrapolations to 1 m?
• I strongly recommend not reporting values that rely on extrapolated depths, and instead presenting results only to the maximum common depth sampled. Spatial variability:
• If the main objective is to assess spatial variability, it is unclear why measurements from different plots, sites, and settings are lumped together.

Response: For each soil parameter (e.g., stocks, Corg content, bulk density, grain size), we have removed the results reporting mean ± SD or SE and included a range for 0-25 cm in the text and provided across plots, sites, and settings in supplementary table Table S2.

• Depth basis revised text. Page 10, Lines 229-230. 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-3 (see site and geomorphology level details in Table S3 and Fig S1).
• We only report measured values in main text and supplements
• Spatial variability: We have addressed this point by reporting spatial variability at different levels through variance partitioning reporting in methods, results and discussion.

Revised manuscript text.

• Methods 2.3. Page 5, Lines 117-123. Our primary carbon stock response was measured-only soil organic carbon stock in the 0–25 cm interval. To ensure comparable depth coverage across sites, primary inference for 0–25 cm stocks was 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.
• Results 3.2 Page 9, Lines 214-221. Sensitivity analyses for 0–100 cm stocks indicated a marginal geomorphology effect for measured-only totals (F (2,2.65) = 8.114, p = 0.074), but no effect after extrapolating all cores to 100 cm (F (2,2.96) = 0.22, p = 0.816). Refusal depth differed strongly among geomorphic settings (Supplement 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 = -77.9 cm; p = 0.002) cores being shallower than estuarine cores (Supplement Table S2). Overall, 25 of 62 cores (40.3%) reached 100 cm depth, and 25 of 62 cores (40.3%) reached ≥80 cm depth.

1. Comment: I would reconsider the discussion around some of the bivariate relationships between carbon stocks and temperature or salinity. The ranges of these variables across sites in this study are extremely small, around 1 deg. C for temperature and less than 3 units of salinity. Given that carbon stocks accumulate over centuries, snap-shot measurements of temperature and salinity are unlikely to explain such large variability in stocks. It is possible that any aren’t relationships are coincidental rather than ecologically meaningful

Response: We have reduced mechanistic claims and interpreted associations cautiously as potentially proxying broader site differences.

Manuscript changes: Discussion 4.3 revised

Revised manuscript text: Page 15; Lines 339-355. Our expanded core-level mixed model found that none of the tested sediment properties (mud fraction, bulk density) or plot-level 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). 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.

1. Comment: The role of seagrass species is not considered at all in the study, yet species composition can strongly influence carbon stocks and sediment properties. I recommend being critical of these results and consider whether these relationships are robust or generalizable beyond your study

Response: We have expanded our discussion to include species/cover/canopy context and note future trait-based modelling.

Manuscript changes: Discussion expanded 4.2 (Page 14, Lines 306–314) + Supplement Table S1.

Revised manuscript text: 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 Corg storage 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., 2025) 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).

1. Comment: L38, L53: The acronym for carbon is used inconsistently (e.g., C_org, C, and carbon). Please select one format and use it consistently throughout the manuscript.

Response: We have standardised notation and terminology across the manuscript.

Manuscript changes: Throughout.

Revised manuscript text:

Old L53. However, typical anoxic conditions and high organic matter inputs in seagrasses that promote high C burial and sequestration can also favour the production and emission of greenhouse gases such as methane (CH₄).

New. However, typical 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₄).

1. Comment: Extrapolating carbon stocks to 1 m depth does not typically lead to underestimation, in most cases it leads to overestimation. Please see Krause et al. (2025) for a discussion of this issue.

Response: We have revised wording to reflect potential under- or over-estimation cited the reference.

Manuscript changes:  Introduction revised

Revised manuscript text:

Old. Page 2; Line 29-31. At the intra-site scale, extrapolating standard 100 cm Corg stocks from shallower cores is a major caveat, potentially leading to 1.5- to 10-fold underestimation (Dahl et al., 2025; Stankovic et al., 2023).

New. Page 2; Line 45-47. 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; Krause et al., 2025; Stankovic et al., 2023).

1. Comment: L55: Regarding CH4 emissions offsetting carbon accumulation. The concept that methane emissions “offset” carbon sequestration is not strictly correct, although it is commonly interpreted this way. Please, see Neubauer (2021) and Neubauer and Megonigal (2015) to understand why this interpretation is incorrect. Ecosystems that accumulate carbon will ultimately produce a negative net radiative balance over sufficiently long-time scales because methane has a much shorter atmospheric lifetime than CO2. Consequently, CH4 emissions do not permanently offset carbon sequestration, this usually only applies to ecosystems that are very short-lived like restored sites.

Response: We have removed the methane “offset” framing from the revised Introduction

Manuscript changes: Introduction revised

Revised manuscript text.

• Old: Page 2, Lines 54-55. Seagrass meadows appear to be negligible to moderate CH4 sources, ranging from 0.3 to 378 μg m−2 h−1 and potentially offsetting ~12–78 times the CO2 equivalents of their carbon accumulation rates (Asplund 55 et al., 2022; Burkholz et al., 2020).
• New: Page 2, Lines 60-63. Seagrass meadows appear to be negligible to moderate CH₄ sources in short-time scale, ranging from 0.3 to 378 µg m⁻² h⁻¹ however ecosystem will ultimately produce a negative net radiative balance over sufficiently long-time scales because CH₄ has a much shorter atmospheric lifetime than CO₂ (Asplund et al., 2022; Burkholz et al., 2020; Neubauer and Megonigal, 2015).

1. Comment: L93 add reference for “(ref)”.

Response: We have added the reference.

Manuscript changes: Methods 2.3 updated (Page 4–5, Lines 107–108). 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.

1. L100: If cores were sampled to 1 m or to depth of refusal, shorter cores should not be extrapolated to 1 m. Instead, carbon stocks should be reported to a maximum common depth across sites. Although global stocks have traditionally been standardized to 1 m depth, this practice is being reconsidered because many deposits, particularly in seagrasses, do not extend to that depth (e.g., Krause et al., 2025). Additionally, sampling depth of refusal itself could represent an informative result. For example, the authors could report whether the depth of refusal varies among geomorphological settings and whether certain environments tend to develop thicker seagrass sediment deposits, which could provide useful information about carbon storage potential across settings.

Response: We have reported refusal depth as an informative result and used measured-only 0–25 cm as primary inference; 0–100 cm measured-only and extrapolated-to-100 cm are sensitivity analyses; refusal depth contrasts are provided in Table S2.

Manuscript changes: Methods 2.3 + Results 3.2 + Discussion revised; refusal depth sheets added in Table S2.

Revised manuscript text: Discussion 4.2. Page 14; Lines 195-302. 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.

1. Comment: L209: The soil organic carbon density results seem somewhat redundant, as C_org content and dry bulk density have already been reported. Consider moving this section out of the main text and instead use the space to add more analytical depth to the results, for example, by including variance partitioning plots (bar charts showing the percentage of variance explained at each spatial level).

Response: We have removed SOC density figure from the main text, added relevant information in the Supplement (Table S2) and focused the main text on variance partitioning across spatial scales.

Manuscript changes: Results revised; Supplement Table S2 updated.

Revised manuscript text: Page 10, Lines 223-227. Variance partitioning indicated strong scale dependence in variability (Supplement 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.

1. Comment: Table 1; The number of cores for stocks in the upper 15 cm is reported as 58, whereas the number of cores for the upper 1 m (including predicted values) is 62. Does this mean that some cores in this study or stock estimates were shorter than 15 cm? Please clarify.:

Response: Yes, three cores were shorter than 15 cm. We have reported depth completeness explicitly by depth interval (e.g., complete-to-15/25/30/100) and presented refusal depth distributions to clarify why n varies.

Manuscript changes: Discussion 4.1 + Table 1 revised; Page 13.

Revised manuscript text: Page 13, Lines 286-289. Table 1. Comparison of seagrass soil carbon stocks (15 cm, 30 cm and 100 cm depth measured only unless indicated) with selected local and global estimates.

1. Comment: L266-267: Table 1 is referenced in a statement indicating that carbon stocks in the upper 30 cm are larger in reef and lagoonal meadows compared with those in estuarine settings. However, Table 1 does not show this result.

Response: We have revised the text to indicate 15 cm soil C_org stocks to corroborate it with Table 1.

Manuscript changes: Discussion 4.1 revised;

Revised manuscript text: Page 14, Lines 317-319. 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

1. Comment: L287-288: Was tidal range measured or assessed for the study sites? Are there differences in tidal range between sites and geomorphic settings, and could these influence carbon stocks or fluxes?

Response: We have clarified that sites share the regional tidal regime; fluxes were measured during low-tide exposure windows; tidal/inundation effects are treated as a limitation; local exposure duration may vary and contribute to variability.

Manuscript changes: Methods 2.1 + Methods 2.7 updated.

Revised manuscript text:

• Methods 2.1. Page 3, Line 86-87. 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.7m, tidal range >2.5 m) and remain fully submerged during high tide.
• Methods 2.7. Page 7, Line 170-172. We acknowledge methodological limitations due to logistical constraints. First, while we targeted high-resolution spatial replication within sites (≤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). As such, we limit data reporting to per hr, not per day or per year.

1. Comment: More information should be provided either in the Methods section or in the Supplementary. In particular, the manuscript should report site and plot coordinates, as well as the number of sediment cores and flux measurements collected at each plot throughout the study period. It is unclear how many times flux samples were collected. Were measurements taken only once per site, or were there repeated sampling events? The timing of the measurements should also be specified (e.g., month, time of day and tidal stage, ebb or flood). Was there a fan in the chambers? Please clarify whether an equilibration period was allowed after placing the dark chamber tops, or whether the initial concentration data were excluded from the regression. The abrupt transition from light to dark can create a short transient period as photosynthesis ceases and the system stabilizes. In addition, please report how many flux estimates were retained and discarded per site, so that readers can assess how many measurements met basic quality-control criteria

Response: We have added supplementary tables with site/plot coordinates, time, and replication counts; clarified chamber mixing, equilibration period, dark-chamber limitation, and QC criterion for flux retention.

Manuscript changes: Supplement Table S1 updated (coordinates) + Methods 2.4 revised (Page 5–6, Lines 124–132).

Revised manuscript text: Methods 2.4. Page 5-6, Lines 130-141. 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-minute equilibration period before T₀ sampling based on our observations from previous study (Iram et al., 2021).  Chambers were rigid PVC with an enclosed headspace height of 30 cm. Lids were sealed using rubber bands and checked for leaks by repeated ambient readings. Measurements were conducted between 10:00 and 17:00 to minimise diel variability. Within each plot, three chambers were deployed at random positions. 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). Before each sample collection, the syringe was pumped 2–3 times to mix the headspace air. 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 r²≥0.8 met, and non-linear estimates were excluded.