Risks to carbon storage from land-use change revealed by peat thickness maps of Peru

Tropical peatlands are among the most carbon-dense ecosystems but land-use change has led to the loss of large peatland areas, associated with substantial greenhouse gas emissions. To design effective conservation and restoration policies, maps of the location and carbon storage of tropical peatlands are vital. This is especially so in countries such as Peru where the distribution of its large, hydrologically intact peatlands is poorly known. Here field and remote sensing data support the model development of peatland extent and thickness for lowland Peruvian Amazonia. We estimate a peatland area of 62,714 km2 (5th and 95th confidence interval percentiles of 58,325 and 67,102 km2, respectively) and carbon stock of 5.4 (2.6–10.6) PgC, a value approaching the entire above-ground carbon stock of Peru but contained within just 5% of its land area. Combining the map of peatland extent with national land-cover data we reveal small but growing areas of deforestation and associated CO2 emissions from peat decomposition due to conversion to mining, urban areas and agriculture. The emissions from peatland areas classified as forest in 2000 represent 1–4% of Peruvian CO2 forest emissions between 2000 and 2016. We suggest that bespoke monitoring, protection and sustainable management of tropical peatlands are required to avoid further degradation and CO2 emissions. Changes in land use threaten the stability of carbon in Peru’s peatlands, which store almost as much carbon as the entirety of the above-ground Peruvian carbon stock but in 5% of the land area, according to maps of the extent and depth of peat.

PMFB, previous mapping was based on relatively small numbers of peat thickness measurements and did not attempt to model and map the spatial variation in peat thickness 2,22 , one of the major sources of uncertainty in the below-ground carbon stock 2 . Rather, the total below-ground carbon stock for the PMFB was estimated by determining the area of different peat-forming vegetation classes (that is, peatland pole forest, palm swamp and open peatland) and multiplying those areas by a mean below-ground carbon stock for each vegetation class. This approach makes several simplifying assumptions 23 : that these three vegetation classes are always underlain by peat, that peat thickness varies more between than within classes and that other land-cover classes (including some wetland ecosystems, such as seasonally flooded forest) do not overlie peat 2,22 . In fact, field observations indicate that these assumptions are no longer valid; in particular, peat thickness varies substantially in space, which includes within single vegetation classes 3,23 . Data-driven maps that more accurately capture the spatial variation in peat thickness and carbon storage, and that cover not just selected study areas but the whole of Peruvian Amazonia, are required to support national and regional peatland conservation planning.
Although Peruvian peatlands are believed to remain largely intact, thus far there has been no quantitative assessment of the GHG emissions that result from land-cover change. Moreover, they face varied and increasing threats, which include agriculture expansion, illegal mining, oil exploration, infrastructure development and the selective felling of the female Mauritia flexuosa palm for commercial purposes 15,[23][24][25][26] . In recognition of these threats, legislation was recently enacted that, for the first time, mandates the explicit protection of peatlands in Peru for climate-change mitigation 27 . To enforce this legislation effectively will depend on a robust mapping of peatland distribution, and on knowledge of the scale and distribution of recent peatland disturbance, none of which is presently available.
Here we present extensive new field observations ( Fig. 1) to test whether previous evidence of a relationship between distance to peatland edge and peat thickness found in other tropical peatlands 3 also applies in Peru. These data are used along with remote sensing imagery to develop the first data-driven models of peatland extent and peat thickness distribution across the whole of lowland Peruvian Amazonia (LPA). We quantified the spatial variation and total peat carbon (PC) stock of these peatlands, and associated uncertainties. Finally, we used these models, along with national data on land-cover change (Geobosques), to map peatland disturbance and estimate the associated CO 2 emissions for the period 2000-2016.
The distribution of peat thickness across LPA is highly variable, with the greatest mean peat thicknesses predicted in the Tigre (232 cm), Marañón (230 cm), Tapiche (234 cm) and Napo (223 cm) basins ( Fig. 2 Table 3 and Supplementary Fig. 5), which gives confidence in our results. We ran two separate peat thickness models: one for the MDD basin and another for the rest of the study area (which contains 97% of the total peatland area). The model that excluded the MDD basin performed better (P < 0.0001, R 2 = 0.66, root mean square error (r.m.s.e.) = 66%; Supplementary  Fig. 5a) than the MDD model (P < 0.0001, R 2 = 0.38, r.m.s.e. = 70%; Supplementary Fig. 5b). We found a significant linear relationship between peat thickness and distance to peatland edge (P < 0.0001, R 2 = 0.13; Supplementary Fig. 6a). This relationship was more significant when the data from the MDD basin were excluded (which gave R 2 = 0.39, P < 0.0001; Supplementary Fig. 6b) and there was no significant relationship between peat thickness and distance to peatland edge within the MDD data (P > 0.1, R 2 = 0.005; Supplementary  Fig. 6c).

CO 2 emissions from land-use change are small but growing
Our analysis of land-use change data shows that a total peatland area of 1,052 km 2 was drained and/or cleared during 2000-2005, which increased to 1,667 km 2 by 2013-2016 (  Tables 4  and 5 for further details). These estimates exclude emissions from areas in which natural peatland vegetation may have been misclassified in 2000 as secondary forest in the land-cover dataset Geobosques (which amounts to 1,353 km 2 ; Supplementary Table  5). These misclassified areas were revealed by visual inspection of a Google map image of the department of Loreto by someone with local expert knowledge (Fig. 3a). For those areas classified as forest in 2000, as accounted for in Peru's 2016 Forest Reference Emission Level report 28 , emissions from peat decomposition represent 0.99-3.72% of the total national CO 2 emissions from LPA forests (that is, from peat decomposition and biomass loss due to gross deforestation; Table 1).

Synthesis and future directions
Our estimate of the total PC stock of 5.38 (2.55-10.58) PgC across LPA is 75% of a recent estimate of the entire above-ground C stock of Peru 29 , and approximately doubles previous estimates of the Peruvian tropical peat stock calculated for the PMFB and the MDD regions only 2,19,22 . Our maps are driven by intensive field sampling which has, for the first time, generated peat thickness data widely across LPA, and which confirms that substantial peatlands extend far beyond the relatively well-studied PMFB. Across the main peat-forming land-cover classes of pole forest, open peatland and palm swamp, above-ground carbon densities (Supplementary  Table 2 23 ) are an order of magnitude lower than the respective PC densities, and total 0.45 PgC (Supplementary Table 2). Summing the above-and below-ground carbon stocks gives a central estimate of 5.83 PgC stored in LPA peatlands.
The quantitative uncertainties around the peatland carbon stock are reduced compared with those of previous studies, even though our study covers an area more than five times greater 2,22 . Future improvements may be gained by collecting field data where they are still lacking, notably the northwest PMFB and parts of the Ucayali (for example, around Pucallpa) and Morona basins. Unlike previous studies 2,22 , our study placed no constraints on which land-cover classes peat can form under, and we predict that around 2% of seasonally flooded forest is underlain by peat. This suggests that the search for peat should not be solely limited to the well-known peat-forming vegetation types of palm swamp, pole forest and open peatland. In addition to land-cover classification maps, we recommend that future fieldwork is informed by examining maps and remote sensing imagery related to hydrology and inundation, such as height above nearest drainage 30 , normalized difference water index 31 and ALOS-PALSAR 32 (where possible with multitemporal images).
Our approach is driven by remote sensing layers with a global coverage and can thus be readily adapted to other regions, provided sufficient field data are available for calibration and validation. Our results call for caution in treating all tropical peatlands as similar, and demonstrate the importance of field data. For example, the distance to peatland edge has been found to correlate with peat  thickness in other regions, such as the Congo basin 3 , and in most of the basins we studied in Peru. However, we found no significant linear relationship between peat thickness and distance to peatland edge for the data in the MDD basin (P > 0.1, R 2 = 0.005; Supplementary Fig. 6c). Householder et al. 19 suggest that this may be because of specific geological conditions in this region: many of the deepest peats in the MDD are often located adjacent to upland (terra firma) terraces, close to the peatland edge. This means that the relationship between peat thickness and distance to peatland edge is more complex in MDD than in other regions. Past research points to geomorphological differences between the northern and southern parts of Peruvian Amazonia 33 ; although floodplains in  northern Amazonia are often wide, rivers in southern Amazonia more often have narrow floodplains confined by terraces. We recommend that new transects should aim to target a range of landscape types (for example, based on elevation maps) and, where possible, should cover the full cross-section of each individual peatland. In spite of this limitation, our random forest regression model for the MDD region performs reasonably well. This study assessed CO 2 emissions that result from peat decomposition due to land-cover change in Peru. Our results suggest that land-cover change in the peatlands of LPA has thus far been restricted to a few hotspot areas, with the largest area of deforestation identified near Pucallpa in the department of Ucayali, an area in which recent ground observations confirm the presence of deforested peatlands 26 (E.N.H.). Access to these peatlands has been facilitated by the development of roads and the increasing demand for land for commercial plantations (for example, oil palm and rice 34,35 ; D. Garcia-Soria, personal communication). Overall, the estimated emissions from peat decomposition remain low in Peru, but our analysis suggests that the annual emissions are increasing. These findings have two implications for the conservation of these ecosystems. First, the low current emissions support the view that the extensive peatland complex of LPA is an emblematic example of hydrologically intact moist tropical forest with a high structural integrity and therefore should be a high conservation priority 23,36,37 . Investment is required to promote the protection and sustainable management of these widespread and extremely carbon-dense ecosystems, before emissions rise over the coming decades 38,39 . Second, the increasing threats and rising emissions from specific land-use transitions in some peatlands mean that it is important to improve the detection of deforestation and secondary vegetation across the full range of peatland forest types, and to make more extensive measurements of GHG emissions associated with specific land-use transitions across the different forest types [7][8][9] .
Taken together, our results indicate a carbon stock within the peatlands of LPA that is three-quarters as large as the entire above-ground carbon stock of Peru 29 but contained within just 5% of its land area. The peatlands also contribute substantial ecosystem and floristic diversity to the Amazon 40,41 . Although our study indicates that these peatlands remain largely intact, they face varied and growing threats 15,35 . Our mapping and carbon stock estimates may be used to support the implementation and enforcement of recent legislation that aims to reduce emissions 27 and should act to encourage national and international investment in monitoring, protection and sustainable management of Peru's peatlands, so that they avoid a similar fate to the heavily degraded peatlands of Southeast Asia 35 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-022-00923-4.

Fieldwork.
Between 2019 and 2021, we collected 445 new GRPs within LPA ( Fig. 1; 294 of these were presented by Honorio Coronado et al. 23 ), which included data on the substrate (that is, peat thickness, where peat is present) and vegetation type (for example, palm swamp). We focused data collection on regions with no existing GRPs and in which peat was believed to be present based on remote sensing imagery (for example, various Landsat 8 (Supplementary Fig. 7) and Sentinel-2 bands), and included the Napo, Putumayo, Tapiche and Tigre river basins ( Fig. 1 and Supplementary Fig. 1), using the only available means of access, that is, via rivers and streams. We also collected new data on peat thickness and carbon concentration from under-sampled peatland ecosystems (for example, peatland pole forest). We made the sampling as spatially representative as possible within the constraints of logistical feasibility, personal safety and accessibility, which are substantial in these remote regions of Peru. The previously published datasets that we incorporated here were also subject to the same constraints.
Where present, peat thickness was measured with an auger or Russian-type peat corer, either along transects perpendicular to the river at intervals of 200-500 m, or at the four corners and centre of the vegetation plots (see below) in which case the value for peat thickness used is the mean of five point measurements. Working along transects that led away from the river and into the peatlands allowed us to sample across wide hydrological and topographic gradients, which included both minerotrophic and ombrotrophic ecosystems. At 91 of these GRPs, we conducted 1, 0.5 or 0.1 ha vegetation plot surveys (collecting floristic data) for the quantitative classification of ecosystem type 23,41 . Additionally, we used 218 previously published GRPs 2,22,43 (24 with floristic data) collected using a similar transect-based sampling strategy in northern Peru and 465 GRPs 19 (148 with floristic data) collected in southern Peru, which amounted to a total of 1,128 GRPs (Fig. 1). Of these, 887 GRPs ( Supplementary Fig. 8) indicated the presence of peat (defined as an organic layer ≥30 cm thick 44 ). Two examples of transects of peat thickness measurements in the Napo basin are shown in Supplementary Fig. 7.
The majority of peat thickness observations do not have corresponding carbon concentration measurements and thus we cannot enforce a precise cutoff in terms of carbon content. However, we visually identified peat and underlying sediments in the field on the basis of their physical properties (for example, colour, structure and texture) and composition (for example, wood, roots and mineral components) 45,46 . At 35 vegetation plots identified by fieldworkers as being on peat, we took sediment samples in the near-basal peat, transition zone and underlying mineral sediment (typically silts or clays) and measured the loss on ignition (LOI) in each to test the visual assessments. The peat, transition zone and mineral samples had mean LOI values of 70, 28 and 13%, respectively (Supplementary Table 6). This gives us confidence that fieldworkers in this region are able to visually identify peat (in this case, soil with an LOI of at least 50%), as there is typically a clear and distinct transition to mineral sediment in Peruvian peatlands.

Map of predicted peatland extent in LPA.
We created a 50-m-resolution map ( Supplementary Fig. 2) of predicted peatland extent in LPA (defined here as the area covered by two of the ecozones recognized by Peru's Ministry of Environment: Ecozone Selva Baja and Ecozone Hidromórfica 47 ). First, we ran a supervised random forest algorithm (200 trees) in Google Earth Engine to predict the distribution of five classes: peat below forest, peat below non-forest (that is, herbaceous vegetation and shrubland), non-peat below forest, non-peat below non-forest and open water. The model was trained and validated (50/50 split of polygons) using peat-thickness measurements and information on the overlying vegetation, and driven using a stack of seven remote-sensing layers, which included two Sentinel-2 indices (normalized difference vegetation index and normalized difference water index 31 ), three ALOS PALSAR-2 bands (HH, HV and HH/HV 32 ), SRTM 30 m digital elevation 48 (Supplementary Table 7) and an extended version of a land-cover classification produced previously 23 (Supplementary Fig. 9; Supplementary Information has further details). The peat below forest and herbaceous vegetation and shrubland categories were amalgamated to form the map of total peatland extent ( Supplementary Fig. 2. We calculated the 5th and 95th confidence interval percentiles for the peatland area using the area and accuracy of each class, applying the method described in Olofsson et al. 49 (equations (9)-(13)), following Draper at al. 2 and recommended by the Global Forest Observations Initiative.

Model of peat thickness distribution.
Testing showed that peat thickness increases with distance to peatland edge (R 2 = 0.13, P < 0.0001; Supplementary Fig. 6), which indicates that the deepest peat is typically found in the centre of a peatland. We thus calculated distance to peatland edge for each model grid using our map of peatland extent. We used the 1,128 peat-thickness measurements as training data, supplemented with points that we assumed to lack peat located along known rivers and urban areas (based on a combination of local knowledge and inspection of Sentinel-2 and Landsat 8 images), which amounted to a final dataset of 1,359 points. The model was run at a 100 m resolution in Google Earth Engine and driven by the stack of remote sensing layers, with two additional layers: distance to peatland edge and height above nearest drainage 30 (Supplementary Table 8).
To robustly test the model performance, we performed a series of validations to account for spatial autocorrelation. Training the model using data only from within the PMFB (n = 717) and testing against data from outside the PMFB in Northern Peru (Napo, Putumayo and upper Amazon basins, n = 155), the model performed relatively well (observed versus predicted peat thickness, P < 0.0001, R 2 = 0.56; Supplementary Fig. 10a). However, the same model (trained using only PMFB data) was unable to predict the variation in peat thickness observed in the MDD basin data (n = 478, P > 0.50, R 2 = 0.00; Supplementary Fig. 10b). For this reason, we decided to run two separate models for the final analysis, one using data only within the MDD basin (n = 477, number of model trees = 100), and the other using all the other data points (n = 867, number of model trees = 50). The model performance was lower for the model with only MDD data (P < 0.0001, R 2 = 0.38, r.m.s.e. = 70%; Supplementary Fig. 5b) than that using all the other data points (observed versus predicted peat thickness, P < 0.0001, R 2 = 0.66, r.m.s.e. = 66%; Supplementary Fig. 5a). We independently validated both models by training each with 80% of the data (randomly selected) and testing with the remaining 20% ( Supplementary Fig. 5c,d).
To account for the uncertainty associated with our estimate of peat thickness distribution, we ran a k-fold analysis as in Rodríguez-Veiga et al. 50 , splitting the data into 1,000 folds and therefore generating 1,000 predictions of peat thickness per pixel. We took the median, 5th and 95th percentiles of the 1,000 predictions to represent our best estimate (Fig. 2a), minimum (Supplementary Fig. 3a) and maximum ( Supplementary Fig. 3b) peat thickness distributions. We subsequently masked the maps of peat thickness distribution using the map of peatland extent ( Supplementary Fig. 2) and thus restricted our model to only regions predicted to contain peat.
Below-ground carbon stock. A dataset of 68 stratigraphic profiles of carbon concentration (%) and dry bulk density (DBD (g cm −3 )) was compiled using data from refs. 2,22,23,43,51 (Supplementary Table 9). This included ten new peat profiles collected as part of this study and described in Honorio Coronado et al. 23 (Supplementary Table 4 of Honorio Coronado et al. 23 ). We calculated PC stock (MgC ha −1 ) from the peat cores by multiplying peat thickness (cm) by DBD and carbon concentration evaluated at regular intervals down the peat profile to the base of the peat. Laboratory conditions varied depending on the study and can be found in the original papers, along with information on protocols. The studies used a variety of standard methodologies to determine the sample carbon concentrations. In line with our definition of peat, we only retained cores in which the peat was ≥30 cm thick, with a mean LOI of ≥50% and those collected using a Russian corer to ensure that DBD measurements were based on a reliable volumetric sample.
We performed a sensitivity analysis to test which of the three components of PC (that, is peat thickness, DBD and carbon concentration) was the most important. Peat thickness was found to be the most important determinant of total PC (P < 0.0001, R 2 = 0.81; Supplementary Fig. 11). We thus used our model of peat thickness distribution to estimate the total PC for each 100 m grid cell and then summed across the entire LPA to produce a total value for the PC stock.
To produce uncertainty bounds for our estimate of the total PC stock, we ran a Monte Carlo analysis which accounted for the uncertainty in each stage of our methodology. We ran 1,000 simulations for PC, constrained using the standard error of the b-estimates from the regression equation (peat thickness versus PC; Supplementary Fig. 11). This was performed twice, once using the 5th and then the 95th percentile distribution of peat thickness calculated previously ( Supplementary  Fig. 3). These 1,000 PC simulations were in turn multiplied by 1,000 simulations of peatland area per grid, constrained by the confidence intervals calculated previously. Finally, the maps of the 5th and 95th percentile of PC stock per grid were summed across LPA to derive the final minimum and maximum uncertainty bounds.
Activity data and emissions from peat decomposition. To estimate changes in forest cover, we used reports of activity data provided by Peru's national monitoring platform, Geobosques 42 Figure 3 shows our predicted peatland map (produced by re-running our model at a 30 m resolution to match the activity dataset) grouping the categories that represent natural vegetation (forest, forest on wetland, wet savannah, water body and non-forest on wetland), secondary vegetation and deforested areas (agriculture, pasture, urban areas, mining areas and bare ground).
Emission factors for organic soils were taken from Chapter 2 of the 2013 Supplement to the 2006 IPCC Guidelines for the National GHG Inventory for Wetlands 6 (IPPC, Intergovernmental Panel on Climate Change). The values range from 7.5 MgC ha −1 yr −1 for secondary vegetation to 9.6 MgC ha −1 yr −1 for deforested peatlands (Supplementary Table 4). These IPCC values are intended to be used for drained peatlands, but peatland disturbance in Peru does not necessarily entail drainage. Nonetheless, undrained secondary forests on peat in Indonesia lose soil carbon (1.4 MgC ha −1 yr −1 (ref. 10 )) at a similar rate to that of shallow-drained plantations (1.5 MgC ha −1 yr −1 (ref. 6 )), and CO 2 emissions in highly degraded undrained peatlands in Peru (for example, degraded Mauritia-dominated palm swamps classified as secondary vegetation, 7.1 MgC ha −1 yr −1 (ref. 8 )) fall within the range of the values of deforested drained peatlands in Indonesia