Quantifying net loss of global mangrove carbon stocks from 20 years of land cover change

Mangrove forests hold some of the highest densities of carbon recorded in any ecosystem, but have experienced widespread deforestation through conversion to aquaculture and agriculture. Alongside deforestation, mangroves have shown simultaneous natural expansion in some parts of the world, and considerable investments have been made into restoration programmes. Here we estimate net changes in the global mangrove carbon stock due to land cover change between 1996 and 2016, using data on mangrove deforestation and forestation, and proportional changes in carbon stock during processes of mangrove loss and gain. The global mangrove carbon stock declined by 158.4 Mt (95% CI = −156.8–525.9 Mt); a reduction of 1.8% of the stock present in 1996. Efforts to conserve and restore mangroves appear to have had some success, and - along with natural forestation - have contributed to relatively low net losses of mangrove carbon stocks over two decades.

Forestation and deforestation was modelled for each pixel within each mangrove patch as a binomial variable, taking into account the probability of false positive errors 2 . None.

Soil carbon density
Published global mangrove soil carbon map 5 .
The mangrove soil carbon density within each polygon was drawn from a normal distribution with mean taken to be the predicted value and standard deviation to be the reported RMSE of prediction 5 .
Recent research has suggested that mangrove soil carbon densities may have been overestimated in the past 6 . We ran a simulation using corrected estimates according to this study ( The meta-analysis relationship used here was general for all blue carbon ecosystems, not specific to mangroves. We conducted a sensitivity analysis using mangrove-specific at values from two case studies 10

Supplementary Methods 1. Sensitivity analyses
We conducted four sensitivity analyses to quantify the impacts of separate methodological decisions on the conclusions of the study. The sensitivity analysis was conducted across 100 bootstrap replicates, for all patches of mangrove gain and loss within the region of Southeast Asia (defined as the Association of Southeast Asian Nations member states, plus Timor-Leste). Only the region of Southeast Asia was analysed because one of the sensitivity analyses required data on the replacement land covers following mangrove deforestation, which are only available with categories comparable to the meta-analysis of rt values in this region 8,12 . While there is now a global mapping of the broad-scale land cover types driving mangrove deforestation 13 , this categorisation does not map onto the commodity-level categories used in the meta-analysis 8 . In addition to the four sensitivity analysis simulations, we conducted a baseline simulation using identical methods to those followed in the main study, to act as a reference point.

Sensitivity Analysis 1
We investigated the sensitivity of the conclusions to systematic errors in the way that mangrove soil organic carbon stocks have been quantified in past research 6 . A recent study by Ouyang and Lee found that the conversion factor used to translate soil loss-on-ignition (LOI) data to organic carbon content (OC) may have been incorrectly defined in past research 6 . The soil organic carbon dataset we used in this study applied a constant LOI/ OC conversion factor of OC = 0.5 × LOI 5 , while the Ouyang and Lee study found that relationship may actually better described as a polynomial line 6 ( Supplementary Fig. 4). It was not feasible to re-analyse the existing literature and repeat a synthesis and mapping exercise of mangrove soil carbon stocks globally. We therefore developed a method to approximate OC, in order to estimate an Ouyang and Lee-corrected soil organic carbon stock density for each patch of mangrove. The following approximation makes a number of significant assumptions, so we proceeded with the previously-published global mangrove soil organic carbon dataset for the purpose of the main study. However, this approximation is conducted as a sensitivity analysis to explore the potential impacts that applying the Ouyang and Lee correction might have on estimation of global mangrove carbon stock loss. As future empirical research incorporates the Ouyang and Lee conversion factor, more data will become available to repeat global mapping and modelling efforts.
Supplementary Fig. 4. Conversion factors between LOI and OC, including those commonly used in previous literature, the relationship used in the global soil mapping analysis we used to provide soil carbon stock data 5 ,and the relationship proposed by Ouyang and Lee 6 .
To apply the Ouyang and Lee correction ideally requires data on the organic carbon content of the soil (OC); information that is not mapped spatially in the global mangrove soil carbon dataset we used here 5 . We therefore developed a method to approximate OC, in order to estimate an Ouyang and Lee-corrected soil organic carbon stock. The carbon stock density (CD) is typically estimated as; With BD being the bulk density of the soil, and CF being the fraction of coarse matter present. We simplify this equation by assuming that CF is equal to zero, and that the bulk density BD can be estimated from the organic carbon content OC, following the empirically-estimated equation 5 ; Combining these equations allows us to solve for OC based a known value of CD; This equation was solved numerically for OC using the uniroot function in the R statistical programming language 14 , for the range of CD values present in the soil carbon stock dataset 5  Once the organic carbon content for the global mangrove soil carbon dataset (OCs) was known for each mangrove polygon, we converted it to the Ouyang and Leecorrected value (OCo). This correction was made by cross-referencing the difference between the OC = 0.5 × LOI line, which was the conversion factor used in the global soil carbon mapping study 5 , and the Ouyang and Lee line at OC = 0.21× LOI 1.12 ( Supplementary Fig. 4). We modelled the difference between the original and Ouyang and Lee values of OC as a generalised linear model using a quasibinomial error structure (Supplementary Fig. 6; Supplementary Table 5). We used this modelled relationship to estimate the soil organic carbon content of each mangrove patch under the Ouyang and Lee correction OCo. We then converted these values back to estimate the soil organic carbon stock density after accounting for the Ouyang and Lee correction (CDo), by reversing the linear regression equation for the line shown on Supplementary Fig. 5. These Ouyang and Lee-corrected estimates of soil organic carbon stock density were then used to re-run the modelling of net changes in mangrove carbon stock. All other parameters were maintained at the values used in the baseline study.
Supplementary Fig. 6. Relationship between the soil organic carbon content estimated using the Ouyang and Lee correction (OCo) and the soil organic carbon content estimated using the OC = 0.5 × LOI line (OCs

Sensitivity Analysis 2 and Sensitivity Analysis 3
We investigated the sensitivity of the study conclusions to different methods of estimating the proportion of mangrove carbon lost following mangrove deforestation (rt). In the main study, we used temporally-varying rt values taken from a systematic meta-analysis 8 . These temporally-varying estimates represent the overall rt values calculated across all types of mangrove conversion. However, the same metaanalysis showed that conversion of mangroves to different replacement land covers can have different impacts on the loss of carbon stocks 8 . Estimates of rt for common types of post-mangrove land cover types, and confidence intervals, are provided in the meta-analysis 8 , but we did not use these values in the study, firstly because they were only available as static mean values, rather than temporally-varying values 8 . Secondly, data on replacement mangrove land cover is not available globally at the required typological resolution 12 .
In Sensitivity Analysis 2 we used static overall mean rt values extracted from the meta-analysis, and in Sensitivity Analysis 3 we used static land cover changespecific rt values. The rt values were randomly generated for each polygon using a uniform distribution within the 95% confidence intervals reported in the metaanalysis 8 . Data on the replacement land covers that followed mangrove deforestation between 2000 and 2012 were extracted from a previously-published study of Southeast Asia 12 . We assumed that all deforestation within each 2 degree grid cell was caused by the most commonly-occurring form of replacement land cover within the grid cell 12 . Grid cells that did not coincide with any recorded replacement land covers were assigned overall mean rt values. The alternative rt values were then used to re-run the modelling of net changes in mangrove carbon stock. All other parameters were maintained at the values used in the baseline study.

Supplementary Analysis 4
We investigated the sensitivity of the study conclusions to an alternative method of estimating the proportion of mangrove carbon accumulated in the ecosystem following mangrove forestation (at). In the study, we used a whole-ecosystem relationship between the time since ecosystem restoration, and the proportion of the reference carbon accumulated 9 . This relationship was estimated from a systematic meta-analysis, but represents a general pattern across blue carbon ecosystems, rather than being mangrove specific. There is no comparable global meta-analysis specific to mangroves, but there are case studies of accumulation of mangrove soil carbon stocks 10 and biomass 11 in foresting mangroves ( Supplementary Fig. 7). Accumulation curves for the proportion of soil carbon 10 and biomass volume 11 present in afforesting mangroves were taken from two case studies from the published literature 11,15 . Biomass volume was used as a proxy of above-and belowground carbon 11,15 . In Sensitivity Analysis 4 we used these alternate at values to investigate the impacts of using mangrove-specific values. It was not possible to incorporate uncertainty in these values. All other parameters were maintained at the values used in the baseline study.
Supplementary Fig. 7. Temporal accumulation of (a) mangrove soil organic carbon stocks and (b) mangrove tree volume.

Sensitivity analysis results and discussion
Sensitivity analysis 1 gave estimates that differed the most from the baseline simulation ( Supplementary Fig. 8), with a 35% lower median estimate. Sensitivity analysis 2 and 3 gave similar results to each other, with medians that were 27% higher than the baseline. Sensitivity analysis 4 gave almost identical results to the baseline estimate ( Supplementary Fig. 8). All of the median estimates given by the four sensitivity analyses were well within the 95% confidence intervals of the baseline estimate, although there is considerable uncertainty surrounding the baseline estimate.
Sensitivity analysis 1 reiterates the argument that past research has substantially over-estimated mangrove carbon stocks, by over-estimating mangrove soil organic carbon content through the use of an inappropriate conversion factor from loss-onignition 6 . For the purpose of the current study, we chose to continue with the previously-published global mangrove soil organic carbon datasets, as our method of estimating OC in order to make a correction factor should be regarded as indicative rather than robust. However, it is likely that the net losses in global mangrove carbon between 1996 and 2016 may be even lower than reported in our study. As future field measurements and meta-analyses begin to incorporate Ouyang and Lee's conversion factor, we may expect new datasets to become available that could be used to re-estimate global changes in mangrove carbon stock, following the framework we developed here.
Sensitivity analysis 2 and 3 gave higher estimates of the loss of mangrove carbon stocks, because they used overall mean rt values rather than taking into account the amount of time since mangrove deforestation. The overall mean rt values synthesised by the meta-analysis were high because most of the studies compared reference and deforested mangrove areas several years after the deforestation events 8 . The rt values used in Sensitivity analyses 2 and 3 are therefore likely to represent the longer-term carbon stock losses that will eventually result from the mangrove deforestation observed between 1996 and 2016. The minor differences between the results of sensitivity analysis 2 and 3 suggest that it is not too critical to account for land-cover dependent rt values, because these values do not differ greatly between replacement land cover types 8 . Sensitivity analysis 4 gave almost identical results to the baseline simulation, indicating that the simulation is not sensitive to the parameterisation of at value. This may partly be due to the broadly similar shape of the at curves used in the baseline simulation and sensitivity analysis 4. Perhaps more significantly, the relatively lower area of forestation means that this parameter has a relatively low weight in impacting the final outcomes for Drt -Fat.
Supplementary Fig. 8. Sensitivity analysis of net loss in Southeast Asian mangrove carbon stock (Drt -Fat ) between 1996 and 2016. Error bars indicate 95% bootstrap confidence intervals.

Supplementary Methods 2. Mangrove gain and loss classification errors
The GMW classification of mangrove forest has an error rate 2 , leading to quantifiable uncertainty over the presence of mangroves at each location in 1996 and 2016. For each pixel of recorded mangrove forestation or deforestation, there is a probability that it is an erroneous or false positive example. For each pixel that is recorded as mangrove in both 1996 and 2016, or non-mangrove in both 1996 and 2016, there is again probability of error -a false negative case of gain or loss. Here we evaluate the potential for these errors to bias estimates of mangrove carbon stock loss and gain.
The false positive error-corrected area of mangrove forestation (Gp) can be estimated following the equation;

Gp = Gr × Mc × Nc
In which Gr is the area of mangrove forestation reported in the GMW dataset, Mc is the probability that a pixel reported in GMW as mangrove was correctly reported as mangrove, and Nc is the probability that a pixel reported as non-mangrove was correctly reported thus.
The false negative error-corrected area of mangrove forestation (Gn) can be estimated as; In which Lr is the area of mangrove deforestation, Sr is the area of mangrove reported in both 1996 and 2016, and Ur is the area of non-mangrove reported, in both 1996 and 2016 in the GMW dataset. The Mm and Nm terms refer to the probabilities that mangrove and non-mangrove pixels reported by GMW were misclassified.
The total error-corrected area of mangrove forestation (Gt) can be summarised as;

Gt = Gp + Gn
Similarly, the false positive error-corrected area of mangrove deforestation (Dp) can be estimated following the equation;

Dp = Lr × Mc × Nc
The false negative error-corrected area of mangrove deforestation (Ln) can be estimated as; And the total error-corrected area of mangrove forestation (Dt) can be summarised as;

Dt = Dp + Dn
The probabilities of mangrove and non-mangrove classification error (Nc, Nm, Mm, and Mc) can be taken from the GMW classification error matrix 2 . As error matrices for all years in GMW are not available, we assume that they all follow the same error rates as 2010, which is the best documented year 2 . The areal extent of Gr , Lr, and Sr can be calculated by cross-referencing the 1996 and 2016 layers of GMW mangrove extent [2][3][4] . However, to our knowledge, the areal extent of non-mangrove classified by GMW (Ur) was not reported in the associated publications and documentation 2-4 , making it challenging to robustly assess the total error-corrected area of mangrove forestation and deforestation. For the purpose of evaluating potential bias introduced by false positive and negative-error correction, we take a conservative estimate of the area of non-mangrove under the GMW study, and assume that the area of Ur is equal to Sr (Supplementary Table 6).
Supplementary The false negative error-corrected estimates of mangrove forestation and deforestation area were very similar, with the ratio of Gn: Dn almost 1:1. This is due to the similarity of the equations used to generate these estimates, as both incorporate the (Sr × Mm × Mc) and (Ur × Nc × Nm) components. The area of mangrove and non-mangrove that does not change between 1996 and 2016 is greatly larger than either the area of gain or loss, so this part of the calculation largely determines the estimate. Conversely, the ratio of Gp: Dp is more than 1:3, indicating that the false positive error rate contributes more to overestimating deforestation than forestation.
In the study, we included uncertainty in mangrove forestation and deforestation due to the false positive error rate in the bootstrap simulation, because of the potential bias that would occur if this were not accounted for. We did not include false negative error corrections due to the lack of Ur data available to estimate these values, and the negligible bias towards either forestation or deforestation area caused by omitting these corrections.