Datasets
We use the Amazon basin (http://worldmap.harvard.edu/data/geonode:amapoly_ivb, accessed 28 January 2021) as our region of study. To determine the grid cells that are contained within Brazil for a subset of analysis, we use the ‘maps’ package in R (v.3.3.0; https://CRAN.R-project.org/package=maps). This is also used in the plotting of country outlines. The main dataset used to determine forest health is from VODCA33, of which we use the Ku-band product. These data are available at 0.25° × 0.25° at a monthly resolution from January 1988 to December 2016. We also use NOAA AVHRR NDVI34. For precipitation data, we use the CHIRPS dataset40 downloaded from Google Earth Engine at a monthly resolution. Finally, to determine land cover types, we used the IGBP MODIS land cover dataset MCD12C1 (ref. 37). All these datasets are at a higher spatial resolution than the VODCA dataset and thus we downscale them to match the lower resolution. Our SST data comes from HadISST49, where we define a North Atlantic region (15–70° W, 5–25° N), for which we take the spatial mean. The mean monthly cycle is then removed to produce anomalies.
For the vegetation datasets that we measure the resilience indicators on (below), we use STL decomposition (seasonal and trend decomposition using Loess)51 using the stl() function in R. This splits time series in each grid cell into an overall trend, a repeating annual cycle (by using the ‘periodic’ option for the seasonal window) and a residual component. We use the residual component in our resilience analysis. The first 3 yr of data had large jumps in VOD which were seen when testing other regions of the world as well as in the Amazon region. Hence, we restrict our analysis to the period January 1991 to December 2016.
To test the robustness of the detrending, we also vary the size of the trend window in the stl() function. The results from these alternatively detrended time series are shown in Supplementary Fig. 4. The results are also robust to varying the window used to calculate the seasonal component rather than using ‘periodic’; at the strictest plausible value of 13, we still see the same increases in AR(1) (Supplementary Fig. 5).
For the AMO index shown in Supplementary Fig. 13, data come from the Kaplan SST dataset and can be downloaded from https://psl.noaa.gov/data/timeseries/AMO/.
Grid cell selection
We use the IGBP MODIS land cover dataset at the resolution described above to determine which grid cells to use in our analysis. The dataset is available at an annual resolution from 2001 to 2018 (but we only use the time series up to 2016 to match the time span of our VOD and NDVI datasets). To focus on changes in forest resilience, we use grid cells where the evergreen BL fraction is ≥80% in 2001. Grid cells are treated as human land-use area if the built-up, croplands or vegetation mosaics fraction is >0%. We remove grid cells that have human land use in them from our forest analysis, regardless of if there is ≥80% BL fraction in the grid cell.
We measure the minimum distance between forested Amazon basin grid cells and human land-use grid cells in 2016 (believing this to be the most cautious and least biased way to measure distance) using the latitude and longitude of each grid point and computing the great-circle distance. We use human land-use grid cells over a larger area than the basin, so that we can determine the closest distance to human land use, regardless of whether this human land use lies within the basin. We also measure the minimum distance from human land use or roads in Brazil, where we have reliable data on state and federal roads (https://datacatalog.worldbank.org/dataset/brazil-road-network-federal-and-state-highways). As in the main text, we reiterate that these minimum distances can be viewed as the maximum distance from human land use as our data will not include roads for the full Amazon basin, or non-federal or non-state roads in Brazil that will have human activity associated with them.
To ensure that the pattern of changes in resilience is not a consequence of more settlements being in the southeast of the region, combined with the gradient of rainfall from northwest to southeast typical of the rainforest, we measure the correlation between MAP and the distances from the urban grid cells, which is very weak (Spearman’s ρ = 0.109, P < 0.001) and as such we are confident that there are separate processes that causes these relationships.
Resilience indicator AR(1)
We measure our resilience indicator on the residual component of the decomposed vegetation time series. We focus on AR(1), which provides a robust indicator for CSD before bifurcation-induced transitions and has been widely used for this purpose23,25,32. We measure it on a sliding window length equal to 5 yr (60 months). The sliding window creates a time series of the AR(1) coefficient in each location. Our results are robust to the sliding window length used, as shown in Supplementary Fig. 6.
From linearization and the analogy to the Ornstein–Uhlenbeck process, it holds approximately that for discrete time steps of width Δt (1 month in the case at hand):
$$mathrm{AR}left( 1 right) = mathrm{e}^{left( { – kappa {Delta}t} right)},$$
where κ is the linear recovery rate. A decreasing recovery rate κ implies that the system’s capability to recover from perturbations is progressively lost, corresponding to diminishing stability or resilience of the attained equilibrium state. From the above equation it is clear that the AR(1) increases with decreasing κ. The point at which stability is lost and the system will undergo a critical transition to shift to a new equilibrium state, corresponds to κ = 0 and AR(1) = 1, respectively.
Measuring AR(1) across the whole time series provides information about the characteristic time scales of the two vegetation datasets we use26. Inverting κ gives the characteristic time scale of the system; for the VOD, we find 1/κ = 1.240 months, whereas for the NDVI, we find 1/κ = 0.838 months when using the mean AR(1) value across the region. This suggests that, in accordance with our interpretation of the two satellite-derived variables, the NDVI is more sensitive to shorter-term vegetation changes such as leaf greenness, while the VOD’s Ku-band is sensitive to longer-term changes such as variability in the thickness of forest stems.
Creation and tendency of AR(1) and variance time series
For analyses where either MAP bands or distance bands are used to create an AR(1) or variance series, we calculate the mean AR(1) or variance value in each month for forested (BL ≥ 80%), non-human land-use Amazon basin grid cells, from which the tendency of this mean series can be calculated. Alternatively, the Kendall τ for each band can be calculated by taking the mean Kendall τ for each individual grid cell that is within the band. Results from the first option are shown in Figs. 4 and 5 and results from the second method in Supplementary Fig. 15 for AR(1).
The tendencies of the CSD indicators are determined in terms of Kendall τ. This is a rank correlation coefficient with one variable taken to be time. Kendall τ values of 1 imply that the time series is always increasing, −1 implies that the time series is always decreasing and 0 indicates that there is no overall trend. Following previous work25,52,53, we test the statistical significance of positive tendencies using a test based on phase surrogates that preserve both the variance and the serial correlations of the time series from which the surrogates are constructed. Specifically, we compute the Fourier transform of each time series for which we want to test the significance of Kendall τ, then randomly permute the phases and finally apply the inverse Fourier transform. Since this preserves the power spectral density, it also preserves the autocorrelation function due to the Wiener–Khinchin theorem. For each time series this procedure is repeated 100,000 times to obtain the surrogates. Kendall τ is computed for each surrogate to obtain the null model distribution (corresponding to the assumption of the same variance and autocorrelation but no underlying trend), from which we calculate a P value by calculating the proportion of surrogates that have a higher Kendall τ value.
Robustness tests
To account for the possibility of human deforestation interfering with the signals we observe (which may not necessarily be detected by the MODIS Land Cover dataset) we also use the Hansen forest loss dataset54 to determine grid cells to remove in an alternate analysis. The original Hansen dataset is at a 0.00025° resolution, 1,000 times higher than the VOD dataset and as such for each VOD grid cell we measure the percentage of Hansen grid cells that show some forest loss over the time period. Note that this dataset does not specify if the observed loss is natural or caused by human deforestation. Excluding any VOD grid cells that contain more than a conservative 5% of lost forest grid cells according to this dataset and running the analysis in the main paper shows similar results. Supplementary Figs. 16–18 are recreations of Figs. 2, 4 and 5, respectively.
The loss of forest resilience observed as increasing AR(1) in both vegetation indices is supported by another indicator of CSD, namely increasing variance28—of both VOD (Supplementary Fig. 19) and NDVI (Supplementary Fig. 20). Variance is more strongly affected by changes in the frequency and amplitude of the forcing of a system and as such results could be biased towards individual events. However, we assessed the precipitation time series for changes in variance and found no relationship with the variance signals of VOD and NDVI (Supplementary Fig. 21). Nevertheless, AR(1) is considered the more robust indicator55. As another test of robustness, we partition the grid cells into those that are in floodplains and those that are not (Supplementary Fig. 3). Floodplain data are part of the NASA Large Scale Biosphere–Atmosphere Experiment (LBA-ECO)56. We also calculate the resilience signals for the C-band product of VOD for comparison (Supplementary Fig. 22). Despite the smaller temporal scale of this product, we still see increases in both AR(1) and variance. To account for a change in the number of satellites used to calculate VOD, for the Ku-band we also recreate the dataset by sampling a single random day per month rather than taking a monthly average, to mimic a constant satellite pass for the whole time period (Supplementary Fig. 23). Although this expectedly affects the absolute values of AR(1) and variance, their relative changes over time remain unaffected. To further test the robustness of our results, we looked for similar signals of resilience change in terms of trends in AR(1) in addition to variance in the precipitation time series used, as a change in the forcing could have an impact on the forest that could mistakenly be interpreted as a vegetation resilience loss. However, as for the variance, there is no clear increase in the AR(1) of precipitation, nor do the spatial patterns of both indicators reveal any relationship between changing precipitation AR(1) and variance and the observed vegetation resilience loss. Hence, we are confident that changes in precipitation forcing are not driving the vegetation AR(1) signals.
Source: Ecology - nature.com
