21st Century drought-related fires counteract the decline of Amazon deforestation carbon emissions | Nature Communications

Sea surface temperature anomalies

Time series of sea surface temperature anomalies (SSTA) for the Pacific and Atlantic Oceans were obtained from the National Oceanic and Atmospheric Administration (NOAA), Earth System Research Laboratory portal. The Atlantic Multidecadal Oscillation index (AMO)20, which describes average anomalies of sea surface temperatures (SST) in the North Atlantic basin, typically over 0–70N (https://www.esrl.noaa.gov/psd/data/timeseries/AMO/), is calculated from the Kaplan SST V2 data set as a departure from 1951 to 1980 time period49. The Multivariate El Niño Index (MEI)21 integrates six variables that describe the status of the coupled ocean-atmosphere system over the tropical Pacific from 30N to 30S. The variables integrated into MEI are: (1) sea-level pressure, (2) zonal and (3) meridional components of the surface wind, (4) sea surface temperature, (5) surface air temperature, and (6) total cloudiness fraction of the sky. All seasonal values are standardized with respect to each season and to the 1950–1993 reference period (https://www.esrl.noaa.gov/psd/enso/mei/table.html). The Pacific Decadal Oscillation index (PDO), calculated for the North Pacific Ocean, poleward of 20N, represents the leading principal component of North Pacific monthly sea surface temperature variability for the 1900–93 period22. NOAA’s National Centers for Environmental Information (NCEI) PDO index, used in this study, is based on NOAA’s extended reconstruction of SSTs version 4 (https://www.ncdc.noaa.gov/teleconnections/pdo/).

Global gridded SSTA map for September 2015 (Fig. 1f) was produced using the monthly SSTA NOAA NCEP (National Centers for Environmental Prediction) EMC (Environmental Modeling Center) CMB (Climate Modeling Branch) Global Reyn_SmithOIv2 product (http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCEP/.EMC/.CMB/.GLOBAL/.Reyn_SmithOIv2/). Anomalies are calculated as the departure from the 1971 to 2000 climatology50. Details of the geographical locations in which the calculations of oceanic indices are based on (Fig. 1f) follow Enfield et al.20 for AMO, Wolter and Timlin51 for MEI and Mantua et al.22 for PDO.

Rainfall and maximun cumulative water deficit

We quantified the intensity and duration of the drought across Amazonia by calculating annual rainfall and maximum cumulative water deficit (MCWD) anomalies (z-scores) for 2005, 2010 and 2015 as the departure from the 2003–2014 mean, normalized by the s.d. (σ), using TRMM 3B43 7A product from 2003 to 2015 (https://disc.gsfc.nasa.gov/datacollection/TRMM_3B43_7.html)23. We have excluded from the long-term mean all the analyzed drought years. We calculated monthly values for both variables based on a average rainfall and MCWD for the whole basin. We also produced pixel-by-pixel mean and standard deviations of monthly precipitation at 0.25° spatial resolution for the Amazon Biome. The cumulative monthly precipitation was estimated in mm month−1 considering a 30-day month for all the data sets. The p-values associated with the z-scores (anomalies) calculated in this study were analized statistically considering the standard normal distribution, where positive or negative anomalies (σ) between: 1.65 ≤ σ < 1.96 are significant at 90% confidence level, 1.96 ≤ σ < 2.58 are significant at 95% confidence level and >2.58 are significant at 99% confidence level.

The pixel-based calculation of the annual maximum cumulative water deficit (MCWD) from 1998 to 2015, correponded to the maximum value of the monthly accumulated water deficit (WD) reached for each pixel within each year17. The WD considers the mean evapotranspiration value obtained by ground measurements in different locations and seasons in Amazonia17. Hence, based on the approximation that a moist tropical canopy transpires ~100 mm month−1, when monthly rainfall (P) is less than this value the forest enters into water deficit. The following rule was applied to calculate the WD for each month (n) on a pixel-by-pixel basis (each located at i column and j line), with evapotranspiration (E), fixed at 100 mm month−1:

$$\begin{array}{l}{\mathrm {If}}\,\:{\mathrm{WD}}_{n – 1{\mathrm{ }}(i,j)} – {\mathrm{ }}E_{(i,j)} + P_{n(i,j)} < \,0;\\ {\mathrm {then}}\,{\mathrm{WD}}_{n(i,j)} = {\mathrm{WD}}_{n – 1 (i,j)} – E_{(i,j)} + P_{n(i,j)};\\ {\mathrm {else}}\, {\mathrm{WD}}_{n(i,j)} = 0\end{array}$$

The MCWD was obtained for each pixel as the negative of the minimum value of WD among all the months in each one of the years. The MCWD is a useful indicator of meteorologically induced water stress without taking into account local soil conditions and plant adaptations, which are poorly understood in Amazonia.

Active fire incidence

The intensity of fire incidence across the Amazon biome within the limits of the Brazilian Amazon was calculated using active fire pixels or hot pixels data, available from the INPE’s Center for Weather Forecasting and Climate studies (CPTEC) fire monitoring system (http://www.inpe.br/queimadas). This data is derived from the MODIS (Moderate Resolution Imaging Spectroradiometer) sensor on board the polar orbiting Aqua satellite, based on its afternoon overpasses24. The monthly and annual number of fires was quantified as active fire density (accumulated number of monthly active fire counts) by summing the daily observations with a nominal 1 km spatial resolution within each grid cell with 0.25° spatial resolution. The thermal anomalies (z-scores) were calculated similarly to the previous data sets.

Burned area mapping by land cover

We produced burned area maps from 2008 to 2012, centered in 2010. This period accounted for one year prior and one year after PPCDAm phase 2 (2009–2011). The analysis using burned area information was designed to capture the constant decrease in deforestation rates (−1773 ± 507 km2 year−1), including a drought year (2010). This 2nd phase consolidated the PPCDAm programme focusing mainly on environmental monitoring and control.

For mapping burn scars we follow the methodology described by Anderson et al9. We used daily surface reflectance products MOD09GA and MOD09GQ as well as 8 day surface reflectance products, MOD09Q1 and MOD09A1, collection 5 from the MODIS data set. The dates of the cloud-free images were selected based on the latest day of the month with nadir view covering the Brazilian Legal Amazon. The burned area maps were based on the linear spectral mixing model (LSMM) applied on three spectral bands9: red (band 1, 620–670 nm) and near-infrared (NIR; band 2, 841–876 nm) reflectance bands from MOD09GQ product and the shortwave infrared (SWIR; band 6, 1628–1652 nm) from the MOD09GA product. From the outputs of the LSMM we used the shade fraction image, which contains the relevant information for mapping burnt areas. We first applied a segmentation procedure using a minimum area threshold of 4 pixels (~25 ha). Subsequently, an unsupervised classification was performed followed by a post classification image edition. The post classification image edition was carried out by a skilled human interpreter using the natural colour composites of the corresponding images for comparison, aiming to minimize omission and commission errors normally produced by any automatic classification algorithm.

Subsequently, burn scar maps were combined with the INPE-TerraClass33 land cover maps (http://www.inpe.br/cra/projetos_pesquisas/dados_terraclass.php) to divide burned scars by land cover classes. For this study we focused on three classes: First, old-growth forest, second, secondary forest and third, pasture. All burn scars occuring in areas mapped as deforested or non-forests in the INPE-TerraClass33 data set were excluded from the analysis.

Our burnt area estimate is likely to be conservative as the use of MODIS data at 250 m spatial resolution can underestimate the area burned by approximately 25% in relation to manually digitized burn scars based on 30 m spatial resolution Landsat images52, 53.

Deforestation rates

The annual cumulative deforested area from 2003 to 2015 was obtained from INPE’s PRODES data set (http://www.obt.inpe.br/prodes/index.php)7. These values are Brazilian Government’s official estimates of annual deforestation in Amazonia. PRODES uses Landsat-like satellites, with 20–30 m spatial resolution, to quantify the complete conversion of old-growth forests into agricultural uses. PRODES only maps polygons with an area greater than 6.25 ha and does not account for deforestation of secondary regenerating forests. Images from Landsat-5 Thematic Mapper sensor are the most used in the programme. To overcome commonly encountered cloud cover problems, however, other sensors such as the CCD from CBERS-2 and CBERS-2B, LISS-3 onboard Resourcesat-1, and UK-DMC2 images can also be utilized.

MOPITT data

The interannual variability of Amazonian biomass burning emissions was analyzed using the record of carbon monoxide (CO) measurements from the MOPITT (Measurements of Pollution in the Troposphere) satellite instrument. The MOPITT instrument incorporates gas-filter correlation radiometers operating in both TIR (thermal-infrared) and NIR (near-infrared) spectral bands54. MOPITT began operations in 2000. Results presented in this manuscript are based on the MOPITT Version 6 Level 3 (gridded) monthly mean data set exploiting both TIR and NIR observations25. This product has been thoroughly validated using in situ CO measurements acquired from aircraft in diverse locations, including four stations in Amazonia25, 55. For each month of the MOPITT mission, a basin-wide CO average was obtained by averaging V6 Level 3 CO total column values over all valid one-degree grid cells within the Amazon Basin boundaries. Monthly mean values were then used to calculate annual means.

CO2 emissions from deforestation and forest fires

Reported data on net CO2 emissions for the Brazilian Amazon biome from 2003 to 2012 was acquired from the report on annual estimates of greenhouse gas emissions (http://sirene.mcti.gov.br/documents/1686653/1706227/Estimativas+2ed.pdf/0abe2683-e0a8-4563-b2cb-4c5cc536c336) developed by Brazilian’s Ministry of Science, Technology and Inovation26. The estimates of Brazilian reference levels of CO2 emissions for the Land Use and Land Cover Change (LULCC) sector follows the IPCC Guidelines for National Greenhouse Gas Inventories46, the Good Practice Guidance and Uncertainty Management in National Greenhouse Gas Inventories56 and the Good Practice Guidance for Land Use, Land Use Change and Forestry57.

To analyze the 13 years contribution of gross CO2 emissions from forest fires, studied here, in relation to deforestation gross emissions estimates26, we first performed an ordinary least square regression between MOPITT CO total column and burned forest gross CO2 emissions. First, we carried out a first-order estimate of burned forest committed gross CO2 emissions (Femission) for the mapped years (2008–2012) based on the following equation 1:

$$F_{{\mathrm{emission}}} = A_{{\mathrm{O}},{\mathrm{S}}} \times B_{{\mathrm{O}},{\mathrm{S}}} \times \alpha \times 3.6667$$

(1)

Where AO,S is the total area burned (km2) for old-growth forests (O) and secondary forests (S) mapped in this study, B is the mean biomass C content of Amazon terra firme forests (BO = 16,000 Mg km−2 and BS = 8000 Mg km−2)58, α is the emission factor for forests affected by fires (α = 0.4)59 and 3.6667 is the conversion factor from C to CO2. For this parameterization we did not considered the biomass from selectively logged forests that may have burned. Aboveground carbon in logged forests are estimated to be 35% lower than in undisturbed forests59. However, the mean biomass value used in our study58 includes selectively logged areas from 2000–2004, minimizing the impact of potential biomass overestimation in the calculation of CO2 emissions. Moreover, the underestimation of burned area by our method, described above, can further couterbalance the overall biomass overestimation effect in the final CO2 emissions estimates from forest fires.

To estimate gross CO2 emissions from forest fires for the whole 2003–2015 period, we then performed a least square regression between MOPITT CO total column (independent variable) and log10 transformed burned forest committed gross CO2 emissions (dependent variable). The resulting equation 2 (n = 5, R2 = 0.95, F = 56.4, p = 0.005), with associated standard error values in parentheses, was used to extrapolate the values for the whole period.

$$\log _{10}F_{{\mathrm{emission}}} = – 5.55\left( { \pm 1.02} \right) + 3.91\left( { \pm 0.52} \right) \times {\rm MOPPITCO}$$

(2)

Subsequently, to estimate gross CO2 emissions from deforestation for the whole 2003–2015 period, we regressed the reported deforestation gross CO2 emissions (Demission) values from 2003–201226 against deforestation rates (D)7, and used the resulting equation 3 (n = 10, R2 = 0.99, F = 2307.1, p ≤ 0.001) to extrapolate the values.

$$D_{{\mathrm{emission}}} = 93.56\left( { \pm 16.79} \right) + 0.05\left( { \pm 0.001} \right) \times D$$

(3)

We expect that the risk of double counting deforestation emissions as forest fire emissions is negligible, because just minimum fractions of the forest area that have burned are later deforested. Results from a previous analysis60 demonstrated that only 2.6% of all burned forests between 1999 and 2008 were deforested by 2010.

Data analysis

To test our hypothesis on the prevalence of the influence of drought over that from deforestation on fire incidence, consequently impeding a direct reduction of Amazonian C emissions from reducing deforestation rates, we first evaluated the spatial and temporal patterns and trends of rainfall, fires and deforestation using the raw values and anomalies. We then evaluated if the previously attested strong correlation between fire and deforestation15 was stable during the three phases of the PPCDAm programme. To corroborate this analysis, we removed the influence of annual deforestation rates on total active fire incidence, by calculating the number of active fires divided by the deforested area, and quantified the significance of the relationship between number of active fires per km2 deforested and time. This analysis was based on the logic that if deforestation is the main driver of fire incidence, the number of active fires per km2 deforested must be kept constant through time, even with the known deforestation decline trend.

Furthermore, we repeated this analysis with the independent MOPITT atmospheric CO data set following the same procedure. We also tested the temporal shifts in the atmospheric CO data normalized by active fires. This analysis followed the logic that if the source of CO2 emission was not related to deforestation but was instead related to management fires from pastures, which has no net impact on the atmospheric CO2 concentration, the amount of CO released per fire event should have decreased through time. This is expected as burning in pastures release much less carbon than burning in old-growth forests or forests being converted (deforested).

To disentangle the potential sources of fires (old-growth forests, secondary forests, pastures and deforestation) and validate previous analyses carried out in this study we analyzed burned area data stratified by land cover classes33 to show changes in the area affected forest fires during droughts and how this information on burned area relates with the atmospheric CO data.

Finally, to demonstrate that areas with high fire activity and forest area burned are not always related to regions with high deforestation activity, we analyzed the temporal patterns of positive active fire anomalies (number of grid cells with positive anomalies) and forest plus secondary forest area burned (total area burned) across the deforestation continuum. Deforestation categories were classified in 20 percentiles corresponding to the area deforested each year, in addition to one class where deforestation was not observed (no deforestation). For each grid cell, we extracted information on the direction of the active fire anomaly and the burned area. For each deforestation class we then counted the number of grid cell with positive fire anomalies and summed the total area burned.

All maps produced in this study were based on publically available data7, 23, 50 using ArcGIS 10.

Data availability

The data that support the findings of this study are all publicly available from their sources. Processed data, products and code produced in this study are available from the corresponding author upon reasonable request.