A global map of mangrove forest soil carbon at 30 m spatial resolution

A harmonized globally representative database (available at: 10.7910/DVN/OCYUIT) was compiled from peer-reviewed literature, grey literature and from contributions of unpublished data from a number of researchers and organizations. Details of database development and a statistical summary of the data are given in the supplemental information available at stacks.iop.org/ERL/13/055002/mmedia.

In order to maximize the utilization of available soil carbon data, we developed a machine learning-based model of organic carbon density (OCD) which models OCD as a function of depth (d), an initial estimate of the 0–200 cm organic carbon stock (OCS) from the global SoilGrids 250 m model (Hengl et al 2017), and a suite of spatially explicit covariate layers (X_p):

$$\begin{equation} \begin{array}{l@{}l@{}l@{}} {\rm OCD}(xyd) = d + {\rm OCS}_{\rm SG} + X_{1}(xy)\\ + X_{2}(xy) + ... X_{p}(xy) \end{array} \end{equation}$$

where OCS_SG is the aggregated organic carbon stock estimated for 0–200 cm depth using global SoilGrids 250 m approach down-sampled from 250 m–30 m resolution, and xyd are the 3D coordinates northing easting and soil depth (measured to center of a horizon). Note here that we model spatial distribution of OCD in three dimensions (soil depth used as a predictor) using all soil horizons layers at different depths, which means that a single statistical model can be used to predict OCD at any arbitrary depth. This 3D approach to modeling OCD reduces the need for making complex assumptions about the downcore trends in OCD, and maximizes the use of collected data.

The derived spatial prediction model is then used to predict OCD at standard depths 0, 30, 100, and 200 cm, so that the organic carbon stock (OCS) can be derived as a cumulative sum of the layers down to the prediction depth for every 30 m pixel identified as having mangrove forest in the year 2000 (Giri et al 2011). Importantly, we found that there is a spatial mismatch between the global mangrove forest distribution (GMFD) of Giri et al (2011) and satellite imagery (figure S2). To best resolve this spatial mismatch, we have adjusted the GMFD by growing all vectors by one pixel (~30 m) and then filtering out any pixel that falls over water by using Landsat NIR band (see SI for more details).

Environmental covariates have been compiled to represent the postulated major controls on OCS in soils generally (McBratney et al 2003) and specifically for mangrove ecosystems (Balke and Friess 2016). Covariates included:

• 1.  
  Vegetation characteristics including percent forest cover (Hansen et al 2013) and Landsat bands 3 (red), 4 (near infrared), 5 (shortwave infrared) and 7 (shortwave infrared) for the year 2000 (Hanson et al 2013) retrieved from http://earthenginepartners.appspot.com/science-2013-global-forest/download_v1.3.html;
• 2.  
  Digital elevation data, which at or near sea-level approximately follows forest canopy height (Simard et al 2006), from the shuttle radar topography mission (SRTM GL1; NASA, 2013) was retrieved from https://lpdaac.usgs.gov, maintained by the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC) at the USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota;
• 3.  
  Long-term averaged (1990–2010) monthly sea surface temperature (SST) averaged into four seasons were generated in Google Earth Engine from NOAA AVHRR Pathfinder Version 5.2 Level 3 Collated data (Casey et al 2010) and downscaled to 30 m resolution using bicubic resampling;
• 4.  
  The M2 tidal elevation amplitude product (FES2012) from a global hydrodynamic tidal model which assimilates altimetry data from multiple platforms was used to represent tidal range at each location. The FES2012 product was produced by Noveltis, Legos and CLS Space Oceanography Division and distributed by Aviso, with support from Cnes (www.aviso.altimetry.fr/);
• 5.  
  Averaged (2003–2011) monthly total suspended matter (TSM) averaged into four seasons estimated from MERIS imagery collected by the European Space Agency's Envisat satellite. Processed and validated TSM data was retrieved from the GlobColour project (http://hermes.acri.fr).
• 6.  
  A mangrove typology map delineating mangroves into estuaries and then either organogenic or mineralogenic based on an analysis of TSM and tidal amplitude data (Zu Ermgassen, unpublished data).

Sea surface temperature, tidal amplitude and TSM are 4 km resolution ocean products and needed to be extrapolated to each pixel containing mangrove forest. Missing values in the sea surface temperature, tidal amplitude and TSM were first filled-in using spline interpolation in SAGA GIS, then downscaled to 30 m resolution using bicubic resampling in GDAL. Including SoilGrids and depth, there were a total of 20 covariates used in building the mangrove OCD model.

The ability of the training points to represent the entire covariate space of the global mangrove domain was assessed by conducting a principal components analysis (PCA) on 15 000 randomly selected points and the 1613 points used in the spatial model. Spatial variables were detrended and centered by subtracting the mean and dividing by the standard deviation (s.d.) before entering into the PCA analysis.

Soil carbon typically varies in highly non-linear ways with depth and across the landscape and as such the ability of standard parametric models to capture this variation is limited (Jardine and Siikamäki 2014, Hengl et al 2017). Here we model the spatial (xyd) distribution of OCD using a machine learning random forest model implemented in the ranger package (Wright and Ziegler 2015) in the R environment for statistical computing (R Core Team 2000). Given the clustered nature of the point data, we have implemented a spatially balanced random forest model design. Model performance was assessed with a 5 fold (Leave-Location-Out) cross-validation procedure where 20% of complete locations were withheld in each model refitting (Gasch et al 2015). The relative importance of using SoilGrids as a covariate was assessed by implementing the cross-validation procedure with and without this variable.

Finally, prediction error of OCS for 0–1 m depth was derived for ±1 s.d. based on the quantile regression approach of Meinshausen (2006) and implemented in R via the ranger package. This procedure is relatively computationally demanding so a random subset of approximately 15 000 points were selected to calculate prediction errors. All modeling was run on ISRIC High Performance Computing servers with 48 cores of 256 GB RAM.

Soil carbon stocks were calculated for the global extent of mangroves for the year 2000 by summing the OCS in each pixel for 1 and 2 m depths. Country level carbon stocks were also calculated for the same depths. Given the fringing nature of mangroves, a global spatial vector data layer was built that allocated the offshore area for each country where mangrove forests can be found. It was derived from the Exclusive Economic Zone for each country. This layer was then dissolved with the onshore areas for each associated country and subsequently used to quantify mangrove OCS tonnage and areal extent.

Potential loss of OCS due to mangrove habitat conversion was calculated between 2000 and 2015 by summing the OCS in mangrove forest pixels which were identified to be deforested. While this analysis cannot distinguish between natural and anthropogenic disturbance, human-driven land use change is believed to be by far the dominant driver of deforestation in mangrove ecosystems (Alongi 2002, Murdiyarso et al 2015). We define deforested using the Global Forest Change dataset (Hansen et al 2013) available online from: http://earthenginepartners.appspot.com/science-2013-global-forest. We chose to use this approach for estimating deforestation instead of using the derived mangrove tree cover loss data produced by Hamilton and Casey (2016) as used by Atwood et al (2017) because the Hamilton and Casey (2016) analysis only considered forested area as area actually covered by trees (i.e. if a 100 ha forest has 80% tree cover then it is counted as 80 ha of forest). In our opinion, this definition mischaracterizes forest area extent. Next, an estimate of the soil carbon emissions associated with land use conversion is needed. The amount of OCS lost can be highly variable and depends on the new land use (Kauffman et al 2014, 2016b, Jones et al 2015) and probably on soil properties. Pendleton et al (2012) used a 25%–100% loss range. Donato et al (2011) used a low estimate of 25% of the OCS in top 30 cm and 75% in top 30 cm + 35% from deeper layers as a high estimate. Expanding on earlier work, Kauffman et al (2017) found that on average 54% of belowground carbon (soil + roots) to 3 m was lost after conversion to shrimp ponds and pastures. Given the limited number of studies comparing soil OCS change with land use change, in this work we adopt the same 25%–100% range as used by Pendleton et al (2012) applied to the first meter of soil. Finally, country level statistics for OCS loss were calculated as described above. All global and country level analyses were performed on the 30 m resolution dataset in Google Earth Engine (Gorelick et al 2017).

Zoom In Zoom Out Reset image size

The random forest model was successful in capturing the major variation in OCD across the mangrove database (figure 1(a)) with an R² of 0.84 and a root mean square error (RMSE) of 6.9 kg m⁻³ compared to the mean OCD value of 29.6 kg m⁻³. Cross-validation results yield an R² of 0.63 and an RMSE of 10.9 kg m⁻³ (figure 1(b)), which is the de-facto mapping accuracy to be expected on the field. Low OCD values were slightly over-predicted and high OCD values were under-predicted (figure 1(b)). The initial OCS prediction from SoilGrids 250 m was the most important variable explaining mangrove OCD. Running the 5 fold cross-validation without and with SoilGrids indicated that this single variable explained improved model performance by nearly 50% (R² increased from 0.42–0.63). Seasonal total suspended matter (TSM), depth of sample, mangrove tree cover, Landsat Red band, sea surface temperature (SST), and tidal range were the next ten most important variables, respectively (figure 1(c)). Quantile regression analysis indicated that the full uncertainty (±1 s.d.) about a mean prediction of carbon stocks to 1 m depth averaged 40.4% of the mean OCS with lowest relative uncertainty in the most carbon-rich mangrove forests (figure S7).

Zoom In Zoom Out Reset image size

Projection of the mangrove OCD model to global mangrove forests revealed the distribution of soil carbon storage in these ecosystems (figure 2). The mean (±1 s.d.) OCS to 1 m depth was 361 ± 136 Mg C ha⁻¹ with a range of 86–729 Mg C ha⁻¹. At the national level (table S1), Bangladesh had the lowest per ha stocks, averaging just 127 Mg C ha⁻¹ followed by China and the nations bordering the Persian Gulf and Red Sea with an average OCS of 214 and 233 Mg C ha⁻¹, respectively. The highest per ha stocks were found in many of the pacific island nations, averaging 505 Mg C ha⁻¹ with much of Southeast Asia ranking well above the global mean.

While the national level comparisons are revealing, by modeling at a 30 m resolution much richer details of potential within forest variation in OCS are seen (figure 2). Mangrove forests dominated by sediment laden fluvial inputs typically have consistently low OCS as seen in the Sundarbans and Madagascar (figures 2(a) and (e)). In non-deltaic mangroves, the model appears to have captured the large zonal variation in OCS that is often observed in field studies (figures 2(b) and (c).

Utilizing the Hanson et al (2013) global deforestation analysis (figure S8), we found that 278049 ha (1.67% of total) of area identified as mangrove habitat in the year 2000 was deforested resulting in the committed emission of 30.4–122 Tg C (111–447 Tg CO₂) from mangrove forest soils due to land use change between 2000 and 2015 (figure 3). The relative rank of nations in terms of loss of mangrove forest area and OCS were often not the same (table S1). Indonesia alone was responsible for 52% of this global loss with Malaysia and Myanmar representing another 25% of the global total loss (figure 3(c)). When visualized as a percent loss from year 2000 stocks, a slightly different pattern emerged (figure 3(d)). Guatemala had the highest percent loss of mangrove OCS (0.9%–6.8%) followed by several southeast Asian nations, but high percent losses were also found in several Caribbean island nations as well as the United States and several west African countries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our new estimate of global mangrove OCS of 6.4 Pg C in the upper meter and 12.6 Pg C to 2 m is largely consistent with past efforts to calculate this value (Donato et al 2011, Jardine and Siikamäki 2014, Sanders et al 2016). However, our estimate is double that of Atwood et al (2017) primarily due to their use of the Hamilton and Casey (2016) estimate of mangrove extent instead of Giri et al (2011). Importantly, by using an environmental covariate model, we have been able to make plausible estimates for regions where no sampling has taken place instead of relying on global mean values (i.e. Atwood et al 2017).

The total amount of soil carbon was similar in our analysis and the most comparable analysis, that of Jardine and Siikamäki (2014), but the spatial distribution of carbon-rich versus carbon-poor mangroves varied substantially. For example, we found much higher OCS levels in West Africa than in East African nations (figure 2 and table S1) but the reverse was found by Jardine and Siikamäki (2014). Large discrepancies were also found for Colombia, Sri Lanka and many of the countries bordering the Red Sea. These differences were most likely driven by lack of data in those regions at the time of the analysis by Jardine and Siikamäki (2014) given that nearly half the data in our database was collected after their study was published. Additionally, our analysis suggested a much larger range in OCS (86–729 Mg C ha⁻¹) compared to 272–703 Mg C ha⁻¹ in the analysis of Jardine and Siikamäki (2014). This difference in range was likely due to the inclusion of more data from sub-tropical and temperate mangroves (figure S3).

The depth trend analysis (figure S6) and random forest variable importance (figure 1(b)) both indicated that depth should be considered in calculation of OCS. For locations that were either stable peat domes or located in estuaries receiving large annual sediment loads, a stable OCD profile distribution would be expected and this was found for many sites (figures S6(a) and (e)). However, where mangroves are growing in a mineral matrix that is receiving only low sediment loads, a decline in OCD may be expected as the carbon inputs from the productive mangrove forest would be concentrated in the surface horizons (figure S6(b)). Still in other cases (figure S6(c)), changes in hydrologic/sediment regimes can lead to irregular depth patterns or even an increase in OCD with depth.

While total area of mangroves was a key determinant of total soil carbon storage, amongst the top 25 mangrove OCS holding nations, there was a nearly even split between nations with smaller area of high soil carbon density forests and those nations with lots of low soil carbon density forests (figure 4). Indonesia was the clear exception to this trend with the largest mangrove holdings which also contain rich carbon stocks resulting in Indonesia alone holding nearly 25% of the world's mangrove OCS.

Compared to terrestrial carbon pools, mangrove forests rank low due to their limited spatial extent. For example in the upper meter of soil, permafrost affected soils are estimated to store 472 ± 27 Pg C (Hugelius et al 2014), tropical forests contain ~188 Pg C, and soils under permanent cropping contain ~150 Pg C (table 1). However, on an equal area basis, mangrove forests on average store more soil carbon than most other ecosystems (table 1). Importantly, our analysis has demonstrated mangrove soil carbon is highly variable and many mangroves actually store fairly modest amounts of carbon in the upper one or two meters of soil. While not the focus of this analysis, it is important to point out that while some mangrove forests store modest levels of OCS in the upper meter of soil, they can have high sequestration rates and conversely carbon-dense mangroves can have low annual sequestration rates (Lovelock et al 2010, MacKenzie et al 2016).

Our spatial modelling framework, in which global predictions were combined with local high resolution images, was successful as the general patterns of carbon variation from SoilGrids 250 m were maintained, while the spatial detail was significantly improved by moving from 250 m–30 m spatial resolution. The initial SoilGrids 250 m OCS prediction (Hengl et al 2017) was based upon machine learning algorithms using 237 covariates that covered the main state factors of soil formation (Jenny 1994)—climate, relief, living organisms (vegetation) and parent material—and was four times more important in predicting mangrove OCD than any of the local covariates. However, the local covariates allowed for a much more refined picture of the spatial variation in OCS within mangrove forests (figure 2) which was not captured in the 250 m resolution SoilGrids 250 m prediction. Covariates related to hydrogeomorphology (TSM and tidal range), as hypothesized, were important predictors of local variation in OCD. Both TSM and tidal range were strongly negatively correlated with OCD suggesting that locations with either high sediment loads or strong tidal flushing do not accumulate large carbon stocks. The mangrove typology ended up being uninformative likely because most of the data that went into this classification was already captured in the model.

Zoom In Zoom Out Reset image size

Covariates related to mangrove biomass (SRTM elevation and Landsat bands) were also important in explaining the local distribution of OCD (figure 1(b)). However, as pointed out by Bukoski et al (2017), it is unclear whether the importance of these vegetation-related data are causal drivers of differences in OCD or they just happen to co-vary in the same way as OCD. To further explore the relationship between forest biomass carbon and soil carbon storage, we extracted aboveground biomass data from Hutchison et al (2014) and compared it to our OCS results (figure 5). While a clear positive trend was found between biomass and soil carbon storage (R² = 0.26), there is clearly a lot of variance especially at lower biomass levels where nearly the full range in OCS can be found.

Depth below the soil surface was an important covariate in modeling OCD distribution (figure 1(b)). This finding was supported by the database depth trend analysis (figure S6) which indicated that a flat depth distribution of OCD would be an incorrect assumption 37%–64% of the time. These findings suggest that the simple scaling performed in figure S5 and in several previous assessments (Bukoski et al 2017, Jardine and Siikamäki 2014, Atwood et al 2017) may not be an accurate estimation of total OCS especially when extrapolating from a surface horizon sample alone.

Our analysis suggests that mangrove soils have lost or are at least committed to losing 30.4–122 Tg C due to the land use conversion that occurred between the years 2000 and 2015 (figure 3). Given that at the global level the rate of mangrove forest lost was consistent over this time period (Hamilton and Casey 2016), we estimated an annual soil carbon emission of 2.0–8.1 Tg C yr⁻¹. This value is significantly lower than previous estimates (Donato et al 2011, Pendleton et al 2012) for two reasons. First, we use remote sensing-based measurements of actual mangrove loss instead of applying a large range of annual conversion rates, which are notoriously variable according to their source (Friess and Webb 2011). Second, we have summed the actual OCS values for each of the pixels where land conversion has taken place (i.e. figure S8) instead of applying a conversion rate to a mean OCS value.

The three nations of Indonesia, Malaysia and Myanmar contributed 77% of global mangrove OCS loss for this time period (figure 3). Despite similar area loss (figure 3(a)), Malaysia lost approximately twice as much soil carbon as Myanmar due to the large differences in carbon density between these two nations (mean OCS = 485 ± 57 versus 245 ± 63 Mg C ha⁻¹, respectively). This comparison highlights the importance of using local OCS values for estimating carbon emissions attributed to mangrove conversion.

Not all land use conversions result in equal loss of OCS. Conversion of mangrove forest to shrimp ponds results in a rapid and near complete loss of carbon in the upper meter of soil (Kauffman et al 2014), as well as losses deeper in the soil profile (Kauffman et al 2017). Conversion to other agricultural uses such as pasture for beef production (Kauffman et al 2016b) and cereal crops (Andreetta et al 2016) also appear to result in large soil carbon emissions. However, mangrove degradation and loss due to over harvesting for fuelwood (Jones et al 2015) or due to natural disturbance (Cahoon et al 2003) likely leads to more moderate emissions as decomposition and erosion exceed new plant carbon inputs.

It is important to note that nearly all available data on OCS loss due to conversion to other land uses come from organogenic mangrove forests. In a mineral-dominated mangrove systems with only a few percent sediment OCC, we would not expect the same level of carbon loss as when peat deposits are drained or removed. In fact, reclamation of deltaic sediments for paddy rice cultivation can lead to increases in OCS (Kalbitz et al 2013), although methane emissions would be expected to increase. Additionally, if mangrove habitat is lost due to deforestation without a change in hydrologic regime, mineral-dominated mangroves can continue to accrete carbon, but at a lower rate than in a system that has additional organic matter inputs from the mangroves themselves (Pérez et al 2017).

While we endeavored to ensure that the model input data was of the highest quality possible, there undoubtedly remain unknown errors in the database which are contributing to model error. Machine learning models are particularly sensitive to outlier values and extrapolation (Murphy 2012). Various research groups use different methods for determining the organic carbon concentration (OCC) of a sample with not all publications reporting whether or not results were corrected for occurrence of inorganic carbon or whether or not roots were excluded before further processing. Bulk density (BD) is a difficult parameter to measure accurately in many soils, and based on our analysis of BD versus OCC (figure S1) some reported data are unlikely to be accurate. We developed a procedure to correct potential BD errors, but a pedotransfer function gives only an approximation of the true value. Given the importance of depth in our models (figure 1(b)), it was unfortunate that so many investigations only report OCS for large depth increments. We suggest that future studies using the common practice of collecting subsamples within larger horizon increments (e.g. Kauffman and Donato 2012, Kauffman et al 2016a) report the specific depth increment of the sample rather than that of the entire core. This additional level of transparency in the data would allow mass-preserving splines (Bishop et al 1999) to be fit through the distinct measurement intervals, resulting in unbiased estimates of OCS.

The largest uncertainty in the input data likely resulted from imperfect information on plot location. Whether, accidental or purposeful (i.e. not wanting to identify exact locations), spatially misplaced data in publications are of limited utility in geospatial applications. All of the covariate data for each mangrove point were selected from spatial layers resulting in the potential for a mismatch between the recorded OCS and the spatial predictors. In this study, we visually inspected all coordinates against Google Earth imagery and our adjusted mangrove domain, and then contacted many authors to seek further information and manually adjusted coordinates when we were confident that the adjustments lead to better spatial location. In the final analysis, 199 soil profiles had to be excluded from analysis because we could not confidently locate these plots within the adjusted mangrove spatial domain.