Regional dry-season climate changes due to three decades of Amazonian deforestation

Region and period of investigation.

The study is performed over the deforested regions of the Brazilian state of Rondônia, which lie between 65° W, 60° W, 13° S and 8° S (Supplementary Fig. 1). The LBA-ECO ND-01 (Large-Scale Biosphere-Atmosphere Experiment in Amazonia—Ecology—Nutrient Dynamics) land cover time series (1984 to 2010)18 and PRODES INPE (National Institute of space research—Projeto de Monitoramento do Desflorestamento na Amazonia Legal) estimates of annual deforestation (since 1988)32 are used to quantify vegetation cover change in this study. Spatial scales of forest clearing in this region have spanned a large range over the past 30 years making it a favourable test bed to study different scale-dependent mechanisms of convective triggering4, 5, 13.

Dry-season precipitation variability and land surface-induced convective processes are important for vegetation adaptation6, 23, 24 and in regulating the timing of the transition seasons19, 20, 21, 22. They may hence affect wet-season arrival19, which in turn may feedback to vegetation functioning. Moreover, in Amazonia, land–atmosphere coupling is generally strongest in the dry season33. Also surface-induced mesoscale circulations are strongest during midday hours7, 8, 10, 13. Accordingly, we focus our analysis on the year-to-year changes in the peak dry-season (June, July, August and September: JJAS) afternoon conditions. GridSat cloudiness is analysed at 14:00 LT and 17:00 LT. TRMM 3B42 3-hourly precipitation data25 are analysed at 17:00 LT because the effect of increase in cloudiness on precipitation is observed in the late afternoon in the TRMM data set and also in the LBA-ECO flux tower network data compilation29. PERSIANN and TRMM 3B43 daily-average products are used to analyse JJAS daily average precipitation.

Cloud cover data and detection algorithm.

The ISCCP GridSat data set14, 34, produced under the NOAA Climate Data Record program by the International Satellite Cloud Climatology Project (ISCCP), has been utilized to carry out the multi-decadal analysis in this study. This data set combines several decades of geostationary data from all over the globe, cross-calibrating between different instruments in time and space (in time between different GOES satellites, and in space between GOES and other geostationary satellites that are a part of the ISCCP project14). We analysed the data at 14:00 LT and 17:00 LT for the months of June, July, August and September between 1983 and 2008 using a standard cloud detection algorithm17 summarized in the Supplementary Information.

Precipitation data.

We primarily use the PERSIANN-CDR global precipitation data set15 for our trend analysis over the study period of three decades. This data set is available as daily total precipitation at 25km spatial resolution. The period used in this study is 1983 to 2015. PERSIANN-CDR employs the PERSIANN algorithm on GridSat-B1 infrared satellite data, and employs the training of the artificial neural network using the National Centers for Environmental Prediction (NCEP) stage IV hourly precipitation data and is then adjusted using the Global Precipitation Climatology Project (GPCP) monthly product. Hence, it is not independent of the GridSat data. We also use TRMM 3B43 global monthly averaged and TRMM 3B42 3-hourly precipitation products25 to identify spatial patterns in the present-day period. Both are available at 25km spatial resolution since November 1997. The period of TRMM data used in this study is 2002 to 2014.

Model set-up and evaluation.

We used the Ocean–Land–Atmosphere Model27 (OLAM) to carry out our numerical experiments. A successor of the Regional Atmospheric Modeling System, OLAM is a variable-resolution global circulation model that uses finite-volume discretization of the full non-hydrostatic, compressible Navier–Stokes equations on a hexagonal grid. Such formulation is essential for the modelling of convective processes in mesoscale phenomena such as studied in this paper. Moreover, the variable resolution of OLAM also facilitates interactions between large-scale atmospheric dynamics with mesoscale processes without introducing errors due to lateral boundary conditions. These features of OLAM make it a suitable tool to study mesoscale phenomena. Land surface processes, cloud microphysics, cumulus convection, radiative transfer and subgrid-scale turbulence are parameterized. For details on parametrizations used see ref. 13.

We carried out 11 different numerical experiments (Supplementary Table 2). We used 3 different land cover maps: simulations prefixed with DEF86 and DEF06 used the land cover observed in Rondônia18 in 1986 and 2006, respectively; simulations prefixed with FORpr represent pristine forest and replace all pasture in the region with evergreen forest. To reduce the number of degrees of freedom that are introduced by an interactive ocean, we ran OLAM as an atmospheric general circulation model with prescribed monthly averaged sea surface temperatures (SSTs)35. We used 3 different SST forcings. Simulations suffixed SST80, SST00 and SSTcl are forced with monthly SSTs averaged between 1980–1989, 2000–2009 and 1971–2010, respectively. The Hadley Center’s SST time series35 used for this study are available as monthly averages that are further averaged over the periods mentioned above to produce 3 averaged annual cycles. These annual cycles are used for the respective experiments. We also considered a simulation that controlled roughness length by running deforestation-like simulations in which the height of pasture vegetation in Rondônia was set to be the same as for evergreen forest (‘dyn’ in the simulation name). A similar deforested simulation was performed to achieve a minimum possible horizontal variation of surface sensible heat flux while maintaining the roughness differences between the two vegetation types. In this experiment the pasture vegetation had the same properties as evergreen forest except vegetation height (‘thrm’ in the simulation name). Finally, a deforested experiment with no regional topographic variations was also performed (‘topo’ in the simulation name) by applying the area-averaged topography within 68° W to 57° W and 16° S to 5° S.

Deforested experiments DEF86SST80 and DEF06SST00 were performed to analyse the role of changing scales of deforestation and SSTs in the observed decadal changes in the hydroclimate of Rondônia. Experiments DEF86SSTcl and DEF06SSTcl were performed to separate the hydroclimatic effects of changing scales of deforestation from the effects of observed decadal variability of SSTs. Comparing DEF06SSTcl-dyn and DEF06SSTcl-thrm with DEF06SSTcl allowed us to isolate the impacts of roughness length. Lastly, the DEF06SSTcl-topo experiment allowed us to isolate the coupled effects of regional topography on the deforestation-induced hydroclimatic changes. The difference fields (DEF-FOR experiments) show the net effect of vegetation cover change from pristine forest to pasture removing spatial features common to both experiments.

Each experiment consists of a time-lagged ensemble of 24 simulations. The ensemble is generated by initializing 24 simulations at intervals of 1 day—intervals starting at 0:00 UTC, 8 June 2004 and ending 0:00 UTC, 1 July 2004. All simulations end at 0:00 UTC, 1 September. All simulation results presented in this study are averaged over the respective 24 ensemble members. We chose to simulate only the months of July and August due to computational constraints. We specifically chose to analyse the month of August because observations show that this month (amongst June, July, August and September) presents the highest resemblance with the JJAS averaged values. This is summarized in Supplementary Table 1, which shows that amongst all months analysed the trend in the August dipole moment is the closest to the trend in JJAS dipole moment, and the spatial correlation between the August averaged and JJAS averaged cloud cover is the largest.

The initial conditions are obtained from the National Centre for Environmental Prediction36 atmospheric fields for June 2004. Soil energy and moisture are initialized using values obtained at the end of 15-year OLAM spinup for June 2004. The initialization is made with these values averaged over a 15° × 15° area around Rondônia. For both simplicity and lack of available data, we prescribed a uniform, constant clay loam soil type around Rondônia. However, soil characteristics can vary as a function of time over deforested areas due to continual grazing and poor management and can also vary on small spatial scales. Despite our simplification, our prescription gave an adequate representation of the observed sensible heat fluxes (Supplementary Table 3).

The grid resolution used in all of the experiments is ~8km over Rondônia, which gradually increases to 64km over northern South America and up to 256km over the whole globe. We chose to simulate at this resolution because both the thermally and dynamically generated vegetation breezes can be captured at this resolution11, 13. The vertical resolution is set to be 200m up to a height of 2km, which increases up to 2km near the model top at 45km. This grid set-up is the same as that used in ref. 13.

In our simulations the pasture and evergreen forest vegetation types differ in mainly four parameters: vegetation height, rooting depth, minimum stomatal resistance and albedo. The values of these parameters for evergreen forest vegetation are: 32m, 5m, 286sm−1 and 0.12, respectively. The corresponding values for pasture vegetation are: 0.32m, 1m, 100sm−1 and 0.18. The modelled vegetation roughness length L, which is the relevant variable for surface to atmosphere momentum fluxes, is dependent on vegetation height H and vegetation total (leaf and stem) area index TAI:

We evaluate our simulations with eddy covariance, meteorology and radiation data collected at two eddy flux tower sites in Rondônia—Fazenda Nossa Senhora (pasture site at 62.36° W, 10.76° S) and reserve Jaru (forest site at 61.93° W, 10.08° S) compiled in the LBA-ECO CD-32 (Carbon Dynamics) Flux Tower Network Data29 between March 1999 and September 2002. LBA data collected after 2002 are now also available (P. Artaxo, personal communication, 6 January 2017) and should be utilized in future studies for model validation. We also use precipitation data from the monthly TRMM 3B43 (ref. 25) and surface radiation fluxes from CERES (Clouds and Earth’s Radiant Energy Systems) surface EBAF (Energy Balanced and Filled) product37. Daytime boundary layer height data reported in ref. 38 are also used to evaluate the simulated daytime development of boundary layer height over pasture and forest in Rondônia. Model evaluation is performed for the numerical experiment DEF06SST00 and is presented in Supplementary Fig. 6, Supplementary Table 3 and Supplementary Information.

Dipole moment vector.

Dipole moment vectors of cloudiness and precipitation occurrence are calculated for each year. Only the pixels within the deforested boundary (neglecting areas above 9.5° S as they form a separate deforested patch) of the corresponding year are used in the calculation; thus, the area used to define the dipoles is different in each year. It is a vector with X and Y components equal to a spatial sum of the individual products of the pixel-wise percentage deviations and their distance from an origin. Percentage deviations from the deforested area mean are used because we are interested in spatial patterns rather than larger-scale trends possibly associated with (multi-) decadal modes of climate variability.

Effect of temporal and spatial biases in satellite data.

As observed in Supplementary Fig. 2a–d, the JJAS averaged albedo and brightness temperature from GridSat have negative and positive trends, respectively. Although the performance of precipitation data sets has been evaluated in different parts of South America26, 39, similar temporal biases in precipitation data may also exist. For example, both TRMM and PERSIANN precipitation data sets have been shown to have systematic biases in estimating the rain rate, rain volume, rain area and rain location of mesoscale convective systems26, 40 in South America with TRMM in general performing better than PERSIANN mostly due to the application of monthly bias correction using ground-based precipitation in TRMM. PERSIANN is also found to overestimate rain area and have larger biases in locating rain centres as compared with TRMM, both compared with ground-based estimates26. Due to this reason PERSIANN shows negative rain rate biases as compared with ground-based data.

For the above reasons we use TRMM estimates for model evaluation. Also we designed our observational analysis metrics based on percentage deviations from area mean to minimize the effect of the above-mentioned temporal biases in the satellite data. Also the cross-validation of our results of spatial redistribution of clouds and precipitation between three, semi-independent, satellite-based products—GridSat, TRMM and PERSIANN—also provides greater confidence in the changes in the spatial patterns detected over the study region. It is observed however that the precipitation percentage deviations, over the deforested area, estimated using PERSIANN precipitation data are usually smaller in magnitude than the corresponding TRMM percentage deviations (Supplementary Fig. 3g and i) indicating the smearing of the precipitation area in PERSIANN mentioned above. It should be noted that this study serves to detect trends in the spatial organization of clouds and precipitation with increasing deforestation and not trends in absolute values of clouds and precipitation.

Statistical robustness of the emerging dipole signal.

Statistical robustness of the dipole in the present day is tested in Supplementary Fig. 5. Probability distributions in Supplementary Fig. 5a–d are generated using 14:00 LT JJAS percentage cloud cover averaged over 1983 to 1990 and 2001 to 2008 and the corresponding land cover maps in 1985 and 2005 obtained from LBA-ECO ND-0118. Both data sets are re-gridded to 2.4km by 2.4km using two-dimensional interpolation for cloud cover and by calculating the fractional deforested area under each coarse grid cell. The joint probability distribution functions (panels e–g) for the periods 1983 to 1992 and 1997 to 2008 are obtained using 2,000 bootstrapped samples from each period and the difference between the samples.

Supplementary Fig. 5e–g shows that there exists a statistically significant systematic spatial redistribution of clouds between the early and present-day time periods. The figures show this systematic spatial distribution in the difference field (Supplementary Fig. 5g), which resembles the southeast–northwest dipole of cloud occurrence as observed in the present-day period (Supplementary Fig. 5f). This analysis corroborates our claim that there is a statistically systematic (as shown by the bootstrapped sample) dipole signal present in the present-day times that is different from the uniform cloudiness observed over highly deforested areas in the early time period.

Supplementary Fig. 5a–d shows, first, that a change in scale of deforestation between the 1980s to 2000s results in a change in the distribution of clouds over the deforested area—change from a unimodal to a strong bimodal distribution. In the early period the clouds preferably occur over moderately deforested grids (with deforested fraction ~0.3), but the majority of the sparsely deforested grids are characterized by downdrafts with low cloudiness. But in the present day, an equal number of grids are populated by high as well as low cloud cover as depicted by the bimodal distribution. Second, in the present-day time period the regions of high and low cloud cover occur at similar levels of deforestation (grids with ~85% of deforestation); that is, despite the amount of deforestation being the same, there is a preference of high cloudiness in some regions and low cloudiness in others—signifying a secondary role of thermal triggering.

Consistent increase in dipole strength in various data sets.

Figure 4 shows that there is a consistent three- to fourfold increase in the dipole strength between the early and present-day time periods irrespective of the data set. The following method is used to calculate the data presented in this figure. Simulated data include the 1,717m relative humidity averaged between 13:00 LT and 18:00 LT and precipitation averaged between 16:00 LT and 20:00 LT. Simulated data are obtained from the month of August from DEF86SST80, DEF06SST00, DEF86SSTcl and DEF06SSTcl. Observed data used to calculate dipole strength includes 14:00 LT and 17:00 LT averaged GridSat cloud and daily PERSIANN precipitation. Observational period 1 includes JJAS 1983 to 1994 (except 1987 and 1988 due to unavailable data) and period 2 includes JJAS 1997 to 2008 (except 1998 due to unavailable data). The standard deviation in simulated data represents variability between 24 ensemble members. The standard deviation in observational data represents inter-annual variability.

Statistical tools.

The ordinary least-squares fit is calculated with the MATLAB in-built function fitlm. The significance of the trend line is tested using the p value derived from the t statistics under the assumption of normal errors. Trend line for the cloud occurrence dipole (Fig. 2) is calculated neglecting years that have at least 50% data missing. The same years are also removed from the trend calculation of the precipitation occurrence dipole (Supplementary Table 1). Correlation analysis is done using the MATLAB in-built function corrcoef. Statistical significance of differences is calculated with a two-tailed t-test at the 1% significance level. The null hypothesis tested in each case is—mean cloudiness or precipitation in each pixel is equal to the deforested area mean cloudiness or precipitation. A pixel-wise 5-day running average is applied to the GridSat cloud occurrence and PERSIANN and TRMM 3B42 precipitation time series before the t-test is performed. Additionally, to access the robustness of the test results, a Z-test is also performed on the sampling distribution of the sample means obtained from 200 bootstrapped samples for each pixel. The bootstrapped sampling distributions were verified to be normal on a quantile–quantile plot. All the results of hypothesis tests reported in Fig. 1 and Supplementary Fig. 3 were further verified to be robust under resampling.

Code availability.

All modelled data and code that support the findings of this study can be accessed from the Princeton University’s DataSpace online repository (http://arks.princeton.edu/ark:/88435/dsp017s75df850).

Data availability.

The data/reanalysis that support the findings of this study are publicly available online at http://www.esrl.noaa.gov/psd (NCEP/NCAR 1), http://dx.doi.org/10.7289/V59P2ZKR (ISCCP GridSat34) https://www.ncdc.noaa.gov/cdr/atmospheric/precipitation-persiann-cdr (PERSIANN-CDR41) with http://dx.doi.org/10.7289/V51V5BWQ, https://pmm.nasa.gov/data-access/data-downloading (TRMM, accessed: 25 November 2014), and http://www.metoffice.gov.uk/hadobs/hadisst (Hadley Center’s SST). The land use maps and flux tower data are distributed by the Oak Ridge National Laboratory at http://daac.ornl.gov and the deforestation rate time series is provided by the INPE PRODES project at http://www.obt.inpe.br/prodes.