Drivers and modelling of blue carbon stock variability

Tidal marshes, mangrove forests, and seagrass meadows are important global carbon (C) sinks, 15 commonly referred to as coastal ‘blue carbon’. However, these ecosystems are rapidly declining with little understanding of what drives the magnitude and variability of C associated with them, making strategic and effective management of blue C stocks challenging. In this study, our aims were threefold: 1) identify ecological, geomorphological, and anthropogenic variables associated with C stock variability in blue C ecosystems; 2) create a predictive model of blue C stocks; and, 3) map regional blue C stock magnitude and variability. We had the 20 unique opportunity of using a high-spatial-density C stock dataset from 96 blue C ecosystems across the state of Victoria, Australia, integrated with spatially explicit environmental data to reach these aims. We used an information theoretic approach to create, average, validate, and select the best general linear mixed effects model for predicting C stocks across the state. Ecological drivers (i.e. ecosystem type or dominant species/ecological vegetation class) best explained variability in C stocks, relative to geomorphological and anthropogenic drivers. 25 Of the geomorphological variables, distance to coast, distance to freshwater, and slope best explained C stock variability. Anthropogenic variables were of least importance. We estimated over 2.31 million Mg C stored in the top 30 cm of sediment in coastal blue C ecosystems in Victoria, 88% of which was contained within four major coastal areas due to the extent of blue C ecosystems (~87% of total blue C ecosystem area). Regionally, these data can inform conservation management, paired with assessment of other ecosystem services, by enabling 30 identification of hotspots for protection and key locations for restoration efforts. Globally, these methods can be applied to identify relationships between environmental drivers and C stocks to produce predictive C stock models at scales relevant for resource management.

The manuscript is well-written and reads smoothly. The introduction provides a good overview of different factors controlling sedimentary OC stocks in vegetated coastal ecosystems that have been identified in literature. The material and methods section gives a clear overview of the study site, the data used and how the different models have been constructed, aided by a figure visualizing the workflow. The results section describes the most important results in a concise way and the discussion section frames the results with respect to existing literature. Overall, this manuscript provides an interesting approach to calculating sedimentary C stocks in blue C ecosystems at a large spatial scale and is well-worth publishing.
Thank you very much.

General comments
My main concern with the current manuscript is that it addresses topsoil C stocks, while it reports on 'blue C stocks' throughout the manuscript, without referring to the topsoil aspect. Emphasizing this aspect is, however, important: it is well known that sampling depth can have a large effect on conclusions drawn on relative differences in C stocks between coastal sediments at different locations or in different ecosystems. For example, depending on ecosystem-specific conditions, C stocks are known to decrease substantially with depth below the surface in certain ecosystems, while in others C stocks remain relatively constant with depth. Therefore, I would invite the authors to stress this aspect more throughout the manuscript: (i) the title would be more informative by including that the study concerns topsoil C stocks and (ii) the discussion should include a section where the implications of only considering topsoil samples is discussed.
We appreciate this point and agree that it is important we make it more clear that we are referring to the top 30 cm of sediments. We have done this throughout the manuscript. First, the title has been changed from Drivers and modeling of blue carbon stock variability to Drivers and modeling of blue carbon stock variability in shallow sediments of southeast Australia. Next, we have altered our wording throughout the entire text and in the section headings to specify "shallow sediment C stocks" in place of "C stocks" or "sediment C stocks". We have clarified this point in the methods also: Sediment C stocks to 30 cm deep (referred to throughout the paper as "shallow sediment C stocks")… We have also added a rationale for our 30 cm measurements in the methods: Though it is common in the literature to sample to 1 m deep in blue C sediments, the sampling protocol used for collecting these data (Ewers Lewis et al., 2018) was designed to maximize spatial coverage of shallow sediment C samples rather than sample entire sediment profiles (which may extend well beyond 1 meter deep). Greater spatial coverage allowed us to test the relationships between a variety of potential drivers and surface sediment C stocks on both fine and broad scales.
We have added a discussion of the implications of only measuring the top 30 cm of sediments to the discussion section: We also aimed to maximize our ability to capture relationships between contemporary drivers and sediment C stocks by utilizing sediment C stock data to only 30 cm deep, a sediment horizon more directly impacted by recent environmental conditions compared to deeper stocks due to age. Based on previously estimated sediment accretion rates in blue C ecosystems in the study region (averaging 2.51 to 2.66 mm year -1 in tidal marshes (Ewers Lewis et al., 2019;Rogers et al., 2006a) and 7.14 mm year -1 in mangroves (Rogers et al., 2006a)), the top 30 cm of sediment represents roughly ~113-120 years of accretion in Victorian tidal marshes and ~42 years of accretion in Victorian mangroves. These time scales suggest sediments depths utilized in this study are more appropriate for assessing the impacts of modern environmental conditions on sediment C stocks compared to meter-deep stocks, which can be thousands of years old (e.g. Ewers Lewis et al., 2019). Using shallow sediment C stocks also allows us to be more confident that the vegetation present now has been there during the time of sediment accretion, unlike deeper sediments that are thousands of years old and for which it is difficult to determine what vegetation, if any, was present at the time of accretion.
The variability in shallow sediment C stocks that could not be explained by our modeling may also be related to the inherent challenges surrounding spatial and temporal matching of driver proxies and sediment C stock measurements; the relationship between shallow sediment C stocks and contemporary environmental settings can be represented more accurately for some variables over others.
Ecosystem type was a relatively powerful predictor of shallow sediment C stock variability in our study and this is likely due, in part, to the direct relationship between vegetation and surface sediments. In most vegetated ecosystems, the majority of underground plant biomass and microbial activity exists within the top 20 cm of soils (Trumbore, 2009). For saltmarsh, it has been demonstrated that the top 30 cm of sediment are directly impacted by current vegetation (Owers et al., 2016). Therefore, using shallow sediment C stock measurements allowed us to take advantage of the direct relationship between vegetation and C stocks to explain variability in surface sediments… … …Modern-day factors influencing vegetation can also have impacts on C stocks deeper than the sediments we measured. The effects of underground biomass on sediment C stocks can extend beyond the top 30 cm, and in fact new C inputs and active C cycling by microbial communities can occur as deep as underground roots extend (Trumbore, 2009). These new C additions (and fluxes) at depth fall outside the general pattern of sediment C decay down-core in vegetated ecosystems (Trumbore, 2009) which has previously allowed for linear or logarithmic regressions to be used to extrapolate 1-m deep C contents from shallow (e.g. 30-50 cm deep) sediment C data (Macreadie et al., 2017a;Serrano et al., 2019). The activity of underground biomass and microbes at depth, when considered over space and time, may account for large C fluxes. The influence of anthropogenic activities, such as land use changes, on these processes via impacts to vegetation may largely go unnoticed based on current methods (Trumbore, 2009), both in this study and in blue C stock assessments on larger scales. We suggest further research to understand the dynamics of active C cycling at sediment depths traditionally considered stable… … …It is important to emphasize here that total sediment depths in blue C ecosystems can vary greatly, and are commonly deeper than 30 cm. Blue C ecosystems can have sediments up to several meters deep (e.g. Lavery et al., 2013;Scott and Greenberg, 1983), suggesting the estimates of C stocks measured here are conservative. In spite of these limitations, surface sediment C stock estimates give us valuable knowledge about the sediment C pool most vulnerable to disturbance and how it may be impacted by environmental drivers.
Although I greatly appreciate that the authors have provided a statement that data is available upon request, I would like to ask the authors to consider publishing the data together with the manuscript, or making it available through an online repository, so references to the data can be made. Open data is becoming increasingly important and has the potential to greatly advance the field of wetland biogeochemistry. Thank you for this suggestion. We have inquired about joining and submitting data to the Coastal Carbon Research Coordination Network Data Clearinghouse (https://serc.si.edu/coastalcarbon), which can host our data in one of the Smithsonian Library's digital repositories and be included in the Coastal Carbon Atlas (https://ccrcn.shinyapps.io/CoastalCarbonAtlas/). The data will be issued a digital object identifier (DOI) that can be included with the published paper for readers to easily access the data.
Additionally, to improve transparency and accessibility of the data both utilized and produced in this study, we have added the following table to the supplements (below) with a reference to it in the methods section: Complete details of data availability for inputs and outputs of our models can be found in supplementary  Table S10.
Any modifications made to these data for producing our models are described in the methods section of this manuscript. Please note the reference to the data produced in this study (R code and model prediction rasters) will be updated in this table. Due to both the large size of the raster files and the intellectual property associated with this work as part of the first author's Ph.D. dissertation, the data will be hosted in an online repository (rather than as a supplement) at the time of publication.

Specific comments L83: allochthonous C can also come from terrestrial or estuarine sources
Thank you we have altered the wording to reflect that although we are referring to C transported into the ecosystem via marine tidal flooding, we recognize that C could have originally come from other sources (including, but not limited to terrestrial or estuarine sources), as follows: In higher elevations tidal flooding is less frequent, providing less opportunity for particles and C to settle out of the water column, resulting in a lower contribution of allochthonous C from marine or other sources compared to lower, more frequently inundated marshes… L149-150: would be good to provide a justification of why only the top 30 cm has been sampled Thank you for this suggestion. We agree and have added a paragraph here in the methods to address this point, as follows: Though it is common in the literature to sample to 1 m deep in blue C sediments, the sampling protocol used for collecting these data (Ewers Lewis et al., 2018) was designed to maximize spatial coverage of shallow sediment C samples rather than sample entire sediment profiles (which may extend well beyond 1 meter deep). Greater spatial coverage allowed us to test the relationships between a variety of potential drivers and surface sediment C stocks on both fine and broad scales.
L154-155: would be good to report on the uncertainty associated with the use of spectroscopic techniques to estimate the 'C contents' of the samples. Any idea about the magnitude of this uncertainty? How were C stocks calculated? Were depth profiles of bulk density collected as well? Please briefly explain this, as this is important for the interpretation of the uncertainty on your results. Thank you, we appreciate this reminder to include a brief explanation here (in addition to referencing the original C stock data paper) to aid in the interpretation of our modelling results. We have added the following text to ensure transparency of this information in the present manuscript: Sediment C stocks to 30 cm deep (referred to throughout the paper as "shallow sediment C stocks") were estimated for 287 sediment cores from 96 blue C ecosystems across Victoria in southeast Australia (Ewers Lewis et al., 2018; Figure 1 (0-2, 14-16, 28-30 cm) within each core. Samples were dried at 60℃ until a consistent weight was achieved, then ground. Dry bulk density (DBD) was calculated as the dry weight divided by the original volume for all samples.

Based on the protocols by Baldock et al. (2013), a combination of diffuse reflectance Fourier transform mid-infrared (MIR) spectroscopy and elemental analysis via oxidative combustion using a LECO Trumac CN analyzer was used to determine organic C contents of all samples. Previous studies
have demonstrated the accuracy of using MIR to estimate organic C stocks of sediments (Baldock et al., 2013;Van De Broek and Govers, 2019;Ewers Lewis et al., 2018). MIR spectra were acquired for all samples, then a subset of 200 representative samples was selected based on a principle components analysis (PCA) of the MIR results utilizing the Kennard-Stone algorithm. Gravimetric contents of organic carbon were measured directly in the laboratory for the 200-sample subset (Baldock et al. 2013). A partial least squares regression (PSLR) was created using a Random Cross Validation Approach (Unscrambler 10.3, CAMO Software AS, Oslo, Norway) and used to build algorithms to predict square root transformed total carbon, total organic carbon, total nitrogen, and inorganic carbon for the entire dataset. The PSLR model was evaluated based on parameters from the chemometric analysis of soil properties (Bellon-Maurel et al., 2010;Bellon-Maurel and McBratney, 2011)

, and the relationship between measured and predicted values was assessed based on slope, offset, correlation coefficient (r), Rsquared, the root mean square error (RMSE), bias, and the standard error (SE) of calibration (SEC) and validation (SEP; see Ewers Lewis et al., 2018 for full details). R-squared values for all square root transformed variables were ≥0.94.
Sediment C stocks were calculated based on Howard et al., 2014. Organic C density (mg C cm -3 ) was calculated by multiplying organic C content (mg C g -1 ) by  ). Linear splines were applied to each core to estimate C density for each 2 cm increment within the 30 cm core, then C densities for each interval (measured and extrapolated) were summed and converted to Mg C ha -1 to estimate total stock down to 30 cm deep for each core location… L237: I would refer to table S4 here; this will help the reader to understand how the models were constructed Thank you, we have added a reference to table S4 here as suggested.

L244: it's not clear from the text how the 'averaged models' were obtained and what these exactly are, please explain this in more detail
Thank you for this suggestion. We want to make sure the generation of the averaged models is very clear, so we have added text to the following portion of the methods section: The dredge products of each global model (i.e. models created from "dredging") were ranked using AICc and the best models (delta AICc <2) were used to produce averaged models (named based on the global model they were generated from, e.g. global model 7 -> dredged and averaged -> averaged model 7). Averaged models were produced using the model.avg function ('MuMIn' package v. 1.42.1;Barton, 2018). The parameter estimates for each averaged model represent the average of that parameter's values from the models in which the variable appeared (from within the subset AICc<2).
To help clarify the rationale for the modelling approach we used, which is better for generating robust predictions when complex predictors are involved but cannot utilize standard methods for generating standard errors, we have added the following text to the beginning of section 2.3: To identify drivers of shallow sediment C stock variability and create the best predictive model of shallow sediment C stocks to 30 cm deep we utilized a multi-step process based on an information theoretic approach and multimodel inference (Figure 3). Traditional approaches have relied on identification of the "best" data-based model; however, information-theoretic approaches allow for more reliable predictions through utilization of multiple models, especially in cases where lower ranked models may be essentially as good as the best (Burnham and Anderson, 2002;Symonds and Moussalli, 2011). Further, information theoretic model selection has been demonstrated to provide significant advantages for explaining phenomena with more complex drivers (Richards et al., 2011). Here, we first looked broadly at our variables of interest by narrowing down to the best models containing all possible variables ("global" models, as explained below) using AICc (Akaike information criteria, corrected for small sample size) to explain the variability observed in the training dataset (70% of total C stock data; Symonds and Moussalli, 2011). From there, we identified which variables within the best global models best explained the observed variability in C stock data in order to remove unnecessary variables from the model equation (through the process of "dredging" and selecting the best subset, explained in detail below). The validity of removing unnecessary variables from the model is supported by the concept of parsimony, which suggests models more complicated than the best model provide little benefit and should be eliminated (Burnham and Anderson, 2002;Richards, 2008). The best subset of models generated from the global models ("dredge products") were selected based on delta AICc<2, which are viewed as essentially interchangeable with the best model (Symonds and Moussalli, 2011). Each subset of best models was used to generate an averaged model, which was tested by generating predictions of C stocks for a reserved (30%) subset of the dataset. The best performing model was used to generate a predictive map of C stocks to 30 cm deep for mapped blue C ecosystems in Victoria.

L298: it is not clear what you mean with 'intercept'
Thank you for this comment. We have clarified the definition of 'intercept' and added some text to clarify the meaning of other model outputs: Parameter estimates from averaged models suggests dominant species/EVC was the most important predictor of shallow sediment C stocks, and was the only variable for which the 95% confidence interval of the estimates did not cross zero (Tables 2 and S7), suggesting a true effect of the variable on observed C stock variability (an estimate that included zero means there is potentially no impact of the variable on C stocks). Specifically, seagrasses P. australis, R. megacarpa, Z. muelleri, and Z. nigricaulis had shallow sediment C stocks significantly different than those of coastal tussock saltmarsh (assigned as the intercept in the model, or baseline dominant species/EVC for which to compare the effect of other dominant species/EVCs on C stocks), while all other tidal marsh EVCs, mangroves, and seagrass L. marina did not. This was confirmed by the ANOVA and Tukey's pairwise comparisons… L333: Would be good to provide a measure of uncertainty on the total calculated C stock, similar to the standard deviations you report on the calculated numbers further down in this paragraph. Thank you for this suggestion. The standard deviation estimates you are referring to in the text were calculated by taking the average of the predicted C values for all individual cells overlapping with each ecosystem's areal extent, then taking the standard deviation. The model predictions are spatially explicit; i.e. a predicted C value is generated for each individual raster cell based on the unique characteristics (i.e. combination of spatial data) of that cell that were included in the averaged model. Therefore, it was not possible to predict a single standard deviation or standard error for the total C stock for the entire state of Victoria in the same way because it was a sum.
Instead, the uncertainty of our predicted C stock estimates is represented in our validation step (e.g. Figure S3). We reserved 30% of our C stock dataset to test the accuracy of the predictions generated by each averaged model (as described in sections 3.2 of the results and displayed in Figure S3). These results showed that our best model explained ~half the observed variability in C stocks (Adj R-sq=0.4881; averaged model 2). We have clarified the errors associated with our model predictions by providing further details on the outputs of our validation step, which utilized the reserved dataset (30% of our original carbon dataset) for assessing the prediction power of each averaged model. Using linear regression, we compared predicted values from each averaged model to actual measured C stock values for reserved dataset. The complete outputs for this analysis have been added to the results section 3.2 ("Model validation"), as follows: We have also been more explicit about the error associated with the averaged model used to generate our state-wide shallow sediment C stock predictions by adding text to the results section on modelled shallow sediment blue C stocks (section 3.3) as follows: We estimated a total of over 2.31 million Mg C stored in the top 30 cm of sediments in the ~68,700 ha of blue C ecosystems across Victoria (Table 4; Figure 5). This estimate is based on predictions from our best averaged model that utilized ecosystem type as the ecological variable (averaged model 7), which explained 46.18% of observed variability in C stock data and had an RSE of 39.29.
Due to the complexity of our multimodel approach, we did not include standard errors in predicted values for the entire region of Victoria due to the combination of using an averaged model and having random effects (there are no compatible R packages, to our knowledge).
L336: Please briefly explain how the standard deviations were calculated. What do they exactly represent? Only the spatial variation within these ecosystems, or also uncertainties related to the model procedures used?
The standard deviations here refer to those related to the mean of predicted C stock values for every raster cell of each ecosystem's mapped areal extent. We have updated the text, as follows, to clarify this: Mean predicted C stocks (±SD) to 30 cm deep for each ecosystem type were 57.96 (±2.90) Mg C ha -1 for tidal marsh, 50.64 (±1.35) Mg C ha -1 to mangroves, and 23.48 (±0.57) Mg C ha -1 for seagrass based on predicted C stock values in all raster cells of each ecosystem's mapped areal extent in Victoria. These C stock values ranged from 23. 33 -291.18, 23.34 -77.81, and 23.33 -73.42 Mg C ha -1 for tidal marsh, mangroves, and seagrass, respectively.

L411: I would suggest changing this title to 'Modelled topsoil blue C stocks'
Thank you, we have changed it to "Modelled shallow sediment blue C stocks" here and in the results to be consistent with the wording in the text body and to be transparent about the depth of sediments considered in the study.   ' We have changed this as suggested, thank you.  ' We have changed this as suggested, thank you.

Figure 4: Please explain in the caption what the error bars represent (standard deviation?) and how they should be interpreted.
Thank you for pointing out that this information was missing from the caption. We have updated the Figure 4 caption to more clearly describe how the figure should be interpreted, including the meaning of the error bars, as follows:  (EVC). Bars are color-coded by ecosystem type: red = tidal marsh, green = mangrove, blue = seagrass. C stocks differed significantly by dominant species/EVC, with higher C stocks in coastal tussock saltmarsh, wet saltmarsh herbland, wet saltmarsh shrubland, mangroves A. marina, and seagrass L. marina (group a) compared to seagrasses P. australis, Z. nigricaulis,and Z. muelleri (group b;ANOVA and Tukey pairwise comparison,F8,284 = 34.80,p < 0.001,. Error bars represent one standard error of the mean.

L40: remove 'our'
We have corrected this, thank you.

L250: variable => variables
We have corrected this, thank you.