Observation-inferred resilience loss of the Amazon rainforest possibly due to internal climate variability

Recent observation-based studies suggest that the Amazon rainforest has lost substantial resilience since 1990, indicating that the forest might undergo a critical transition in the near future due to global warming and deforestation. The idea is to use trends in a lag-1 auto-correlation of leaf density as an early-warning signal of an imminent critical threshold for rainforest dieback. Here we test whether the observed change in auto-correlations could arise from internal variability using historical and control simulations of nine sixth-generation Earth system model ensembles (Phase 6 of the Coupled Model Intercomparison Project, CMIP6). We quantify trends in the leaf area index auto-correlation from both models and satellite-observed vegetation optical depth from 1990 to 2017. Four models reproduce the observed trend with at least one historical realization whereby the observations lie at the upper limit of model variability. Three out of these four models exhibit similar behavior in control runs, suggesting that historical forcing is not necessary for simulating the observed trends. Furthermore, we do not observe a critical transition in any future runs under the strongest greenhouse gas emission scenario (SSP5-8.5) until 2100 in the four models that best reproduce the past observed trends. Hence, the currently observed trends could be caused simply by internal variability and, unless the data records are extended, have limited applicability as an early-warning signal. Our results suggest that the current rapid decline in the Amazon rainforest coverage is not foremost caused by global warming.

1 Introduction

The resilience level of the Amazon rainforest to external stresses, such as global warming and deforestation, is subject to ongoing debate (Feldpausch et al., 2016; Boers et al., 2017; Boulton et al., 2022). Since the 1970s the Amazon rainforest has lost about 20 % of its coverage (Simmons et al., 2019), and the net carbon uptake of the formerly persistent carbon sink has been declining over the last 4 decades due to intensification of the dry season and elevated deforestation (Gatti et al., 2021). The rapid development has raised warnings that the Amazon rainforest is approaching a critical threshold, beyond which irreversible damage is unavoidable (Brando et al., 2014; Boers et al., 2017; Boulton et al., 2022; Parry et al., 2022; Doughty et al., 2023).

Model studies, however, differ widely in their results, and the inter-model spread of vegetation responses in future projections remains high with moderate forest resilience this century but a higher risk of sporadic rainforest loss past 2100 (Huntingford et al., 2013; Boulton et al., 2017; Chai et al., 2021). Notably, Parry et al. (2022) found localized rainforest dieback using an abrupt-shift-detection algorithm in five out of seven CMIP6 models they investigated in simulations wherein CO₂ is increased by 1 % yr⁻¹ until it quadrupled after 140 years. All in all, though, climate models do not predict an imminent and complete collapse of the Amazon rainforest. However, it has been reported that climate models underestimate vegetation-related feedbacks (Richardson et al., 2013; Green et al., 2017; Forkel et al., 2019).

Observation-based studies of the recent historical record convey a more alarming picture. Tao et al. (2022) reported the capacity of undamaged rainforests to withstand future droughts to be limited, especially in the Amazon. Analyzing remotely sensed vegetation data, Boulton et al. (2022) supported this idea and found that more than three-quarters of the Amazon rainforest has been losing resilience since the 2000s, especially in regions of less rainfall and in proximity to regions of human activity. They present evidence for an imminent tipping point of the rainforest in the near future.

Such a tipping point may be initiated by a major tree loss from fires, deforestation, or climate change (Cox et al., 2008; Brando et al., 2014). Land cover transitions, such as forest to crops or forest to pasture, decrease the net surface radiation and latent heat flux while increasing the sensible heat flux, resulting in warming of the land surface (Silvério et al., 2015). Reducing the vegetation density by deforestation is associated with enhanced precipitation run-off and reduced evapotranspiration. Hence, both deforestation and forest degradation by droughts weaken the moisture transport by recycling, which is mainly directed westwards over the Amazonian basin along the prevalent wind direction (Salati et al., 1979). This causes reduced precipitation downwind and degraded forest health in a positive feedback loop.

Negative feedbacks and stabilizing effects may also exist. For instance, vegetation responds positively to increasing levels of CO₂ (Kolby Smith et al., 2016), provided sufficient water and nutrients are available, something which can be observed to happen in most parts of the world, including the Amazon basin (Zhu et al., 2016). Another possible mechanism could be convective clouds that actively shift precipitation from wet to dry regions; the temperature gradient, arising from evaporative cooling in wet regions while warming dry regions activates a low-level breeze that transports moisture to the dry areas (Hohenegger and Stevens, 2018). This, by extension, would be a negative feedback in the Amazon in that the atmosphere acts to moisten dry regions. Since the Amazon rainforest has existed for at least thousands (Malhi et al., 2004) or even millions of years (Maslin et al., 2005), the rainforest must have been dominated by negative feedback in the past.

Tipping points are typically accompanied by a regime shift from a stable state, where negative feedback mechanisms dominate, to a marginally stable state with transition to a net positive feedback parameter. Several statistical metrics, known as early-warning signals, have been proposed to predict a regime shift (Scheffer et al., 2009; Lenton et al., 2012). Most commonly, they quantify the recovery rate of the system to small perturbations. The resilience of a system is subsequently defined as the ability to recover from those disturbances. As the stability of the system decreases, it recovers slower when stochastically forced. This phenomenon, known as critical slowing down (CSD), can be detected by an increase in the lag-1 auto-correlation (AR(1)) of a time series representing the dynamics of the system. Increasing AR(1) has been widely used as an early-warning signal for Earth's dynamical systems such as the western Greenland Ice Sheet (Boers and Rypdal, 2021) or the Amazonian rainforest (Boulton et al., 2022). However, the increase in AR(1) has been shown to occur, likewise, for other physical reasons (Verbesselt et al., 2016), as well as to not occur prior to critical transitions; for instance, when the rate of forcing is higher than the intrinsic response timescale for CSD (Boulton et al., 2013).

Here, we compare the observational record to nine large CMIP6 historical model ensembles and control simulations, quantifying model forest resilience between 1990 and 2014 using the same method as in Boulton et al. (2022). By analyzing large model ensembles, we can test whether internal variability in auto-correlation could be the source of the observation-inferred resilience loss.

2 Methods

2.1 Data

We use the Amazon basin as our region of study, taken to be the domain defined as by RAISG (Amazon Network of Georeferenced Socio-Environmental Information; accessed in March 2023). The observational data are provided by the Vegetation Optical Depth Climate Archive (VODCA) (Moesinger et al., 2020), which is available in a 0.25°×0.25° resolution in daily frequency for the period July 1987 to June 2017. The passive or active satellite observations capture the attenuation of microwave radiation by vegetation, which is known as the vegetation optical depth (VOD). The attenuation depends on various factors like the density, type, and water content of the vegetation and the wavelength range of the sensor (Owe et al., 2008). Shorter wavelengths are more sensitive to the upper-leaf canopy than longer wavelengths, since they experience higher attenuation by vegetation. We chose the lowest-wavelength product available (Ku band; ∼ 19 GHz) for the period January 1990 to December 2017, following the study by Boulton et al. (2022). The monthly means of the VODCA product are interpolated to 1°×1° to compare them better to models, although this did not substantially affect the observed trend (see Fig. A1).

Figure 1 Anomaly of the spatially averaged AR(1) series for observations and MPI-ESM1-2-LR piControl (a) and historical (b) ensemble. The anomaly is computed by subtracting the temporal mean from the time series and is plotted at the end of the 5-year sliding window. The ensemble member with the largest and smallest τ_K values are highlighted. The color gradient of all ensemble runs corresponds to increasing τ_K.

To assess internal variability, CMIP6 model ensembles with at least seven historical runs available and interactive leaf area index (LAI) are included in the study (Table 1). The non-dimensional LAI is defined as the total area of leaves per unit surface area. Although LAI and VOD are not identical variables but physically closely related, changes in both variables are strongly correlated, such that they can be considered good proxies for forest health and resilience (Moesinger et al., 2020).

The model LAI output is evaluated in the period January 1990 to December 2014, as the historical experiments are only available until this point in time. Additionally, we use Shared Socioeconomic Pathway 5 (SSP5-8.5) and pre-industrial control simulations of 500 or 1000 years in length that we cut into windows of 25 years that correspond to the length of the historical period. In that way, we create a control ensemble of 20 or 40 members. Control simulations have the advantage that the internal variability in the model can be directly assessed since no external forcing is present.

The models' land surface components (Table 1) simulate exchanges of energy, water, and carbon between the land surface and the atmosphere, as well as biogeochemical processes like photosynthesis. The LAI is computed dynamically in all models investigated. The ecosystem processes are mostly run on a daily time step; however, CanESM5 and IPSL-CM6A run a water and energy budget, as well as photosynthesis computations, on a sub-daily scale (in accordance with the atmospheric component). Only the two MPI-ESM1-2 models and EC-Earth3-Veg have simulated natural vegetation distribution, while all other models investigated prescribe the plant functional type (PFT) fractions e.g., using the LUH2 v2h data set (Hurtt et al., 2020). However, the number of implemented PFTs differs widely from 4 (CanESM5) to 22 (CESM2). The nitrogen cycle is explicitly resolved in ACCESS-ESM1.5, CESM2, MIROC-ES2L, EC-EARTH-Veg3, and MPI-ESM1-2, while wildfires that are deliberately modeled in CESM2 and the MPI-ESM1-2 family occur randomly in EC-EARTH-Veg3 and are not considered in the other models (references see Table 1). Song et al. (2021) report a generally good global representation of LAI in the global climate models (GCMs) in question. In the Amazon, the multi-year (1982–2014) and yearly average is slightly underestimated by models, just as the inter-annual variability is also underestimated. Long-term LAI trends are generally captured but can vary in bias compared to observations in forested regions.

Last, we determine the land use and land cover change (LUCC) using the land use harmonization data set LUH2 v2h and v2f (Hurtt et al., 2020) that is available as annual values on a 0.25°×0.25° spatial resolution and utilized to force the land components of CMIP6 models.

Figure 2 Maps and histograms of the τ_K values of selected members from the piControl (a, c) and historical (b, d) MPI-ESM1-2-LR ensemble. The observed τ_K distribution is plotted in grey for comparison. We selected the ensemble members that have the highest and lowest τ_K values of the spatially averaged AR(1) series to present the full range of the model variability.

2.2 Resilience indicator AR(1)

The random variability in a signal contains information about the recovery rate from stochastic perturbations, such that we can separate trend, seasonality, and residual of the signal using seasonal trend decomposition (STL) by loess (Cleveland et al., 1990). Assuming the seasonality is constant in time, we choose the STL input parameters trend = 19, season = 13, and period = 12 months, corresponding to choices made in previous work (Boulton et al., 2022). Nevertheless, slightly altering the STL parameters has no relevant influence on the presented results. The residual component of each grid cell can then be used to quantify the short-term responses of the forest by calculating the AR(1) on a sliding window of 5 years. We find the AR(1) time series using an ordinary least squares fitting method for the autoregressive model (Eq. 1), where ϵ_t represents the white noise of the model with zero mean and constant variance ${\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}}$, X_t is the time series in each grid point, and φ is the autoregressive coefficient, also denoted AR(1).

$$\begin{array}{}\text{(1)}& X_{\mathrm{t}}=\mathit{\phi }{X}_{\mathrm{t}-\mathrm{1}}+{\mathit{ϵ}}_{\mathrm{t}}\end{array}$$

We then quantify the trend of this AR(1) time series with Kendall's rank correlation τ_K that measures how well two data vectors agree on their ranks (Kendall, 1938). Choosing one vector to be time, a Kendall's τ_K=1 indicates a strictly increasing AR(1) trend, ${\mathit{\tau }}_{\mathrm{K}}=-\mathrm{1}$ a strictly decreasing, and τ_K=0 no trend. The significance p of this statistic is computed by randomly generating phase surrogates of the time series' Fourier transform under constant variance and serial correlation (Dakos et al., 2008).

For comparison between observations and models, we calculate both the τ_K of the spatially averaged AR(1) series and for individual grid cells. The spatial distributions of τ_K are tested for the similarity in the frequency between the observations and the model using the non-parametric, two-sample Kolmogorov–Smirnov test (Berger and Zhou, 2014). Its test statistic is given by the maximum difference between two cumulative distribution functions and is computed under the null hypothesis that both samples are drawn from populations that have an arbitrary, yet identical, underlying distribution. The higher the p value of this test, the less likely it is that the samples are drawn from different underlying distributions.

Ziehn et al. (2020)Boucher et al. (2020)Mauritsen et al. (2019)Hajima et al. (2020)Swart et al. (2019)Mauritsen et al. (2019)Volodin and Gritsun (2018)Danabasoglu et al. (2020)Döscher et al. (2022)

Table 1 CMIP6 model ensembles used within this study.

3 Results

The observational record of the auto-correlation of VOD from the Amazon basin exhibits variations and trends (Fig. 1), but it is impossible based solely on a single-data record to figure out whether such variations are caused by external forcing or whether it is simply an expression of internal variability. A commonly used method to detect forced changes and events in observational records that contain internal and natural variability is to compare them in various ways to multiple climate model runs with different codes and/or starting from different initial conditions (e.g., Hasselmann, 1997; Otto, 2023). In particular, it has become common that global climate models are run multiple times with the same historical boundary conditions but start from different initial conditions in order to explore their internal variability (Kay et al., 2015; Maher et al., 2021; Table 1). With such large ensembles it is possible to ask whether the observed trend is within the range of variability exhibited by the model and by extension, in the present case, whether an increasing trend in auto-correlation constitutes a skillful early-warning signal.

Figure 3 Cumulative distribution functions of τ_K for observations and members of MPI-ESM1-2-LR. The KS statistic is drawn as a black error bar. Members that are not significantly different from the observational distribution are marked in red. Member 28 shows the best agreement with KS =0.08 and p=0.22.

As an example, we display 30 simulations of the historical experiment of one CMIP6 model and MPI-ESM1-2-LR, together with the observations (Fig. 1; right panel). We see that the observed trend in the spatially averaged AR(1) of τ_K = 0.62 is within the range of trends exhibited by the model ensemble members (−0.76 to 0.72) when calculated over the same period. Therefore, in terms of this model's behavior, the observed trend is within the range of variability. To investigate whether global warming or land use change is affecting the trends, we can also inspect equally long-duration chunks from the same model's pre-industrial control simulation (Fig. 1; left panel). Here we also find trends that encapsulate the observed trend (−0.68 to 0.79). Even though the shorter-term deviations appear slightly muted in this case compared to the historical ensemble, the overall trend range in the control and historical experiment does only differ marginally. Moreover, the largest AR(1) trend of this model is found among the control runs, suggesting that the forced response of the Amazon rainforest is not needed to generate an increase in AR(1) of a similar magnitude to that observed.

The modeled spatial distributions of the trends in AR(1) further support the idea that the origin is internal variability. We can inspect maps and frequency distributions from the two historical ensemble members and the two chunks from the pre-industrial control in Fig. 2 that exhibited the largest trends (Fig. 1). All historical ensemble members are shown in Fig. A2. We see that positive trends can occur in any part of the Amazon basin and are not preferentially occurring in the southern parts where most of the land use changes happened (Fig. A6). There are, however, cases wherein this pattern occurs (e.g., members 4, 13, 15, 26, and 30). Also noteworthy is that the model does display trend patterns with spatial scales that are substantially larger than the model resolution such that the underlying causes must be a simulated feature of the model, e.g., large-scale weather events.

Figure 4 Model's mean τ_K distributions for piControl and historical runs. The histograms are fitted with a β distribution on the interval $[-\mathrm{1},\mathrm{1}]$. The observation's percentile within the respective distribution is denoted by p_cml. Historical forcing makes the occurrence of the observations more likely only in MIROC-ES2L and ACCESS-ESM1-5, while CanESM5 and MPI-ESM1-2 simulate a higher probability in control runs.

Figure 5 Historical and SSP5-8.5 scenario of MPI-ESM1-2-LR LAI and the corresponding spatially averaged AR(1) trend compared to the observational VOD record. Panel (a) displays the STL trend component of LAI and VOD respectively, and panel (b) shows the AR(1) series of the STL residual plotted at the end of the 5-year sliding window. Ensemble member 28 is the best-agreeing member with observations in the historical period, while member 10 has the least agreement. Note that the dual vertical axes used in the upper panel are scaled to have the same relative range.

Not all models are equally fit for the purpose of simulating the Amazon rainforest dynamics in accordance with observations. To investigate this, we apply the Kolmogorov–Smirnov (KS) test (Sect. 2.2) on each ensemble member to test whether it could have been drawn from the same underlying distribution as the observations. The process is illustrated in Fig. 3. For the MPI-ESM1-2-LR model, the test identifies two ensemble members (5 and 28) that are statistically indistinguishable from the observed trend distribution. This implies that the observations lie within the range of model variability but on the edge of what the model is able to reproduce. We carry out the same procedure for all models and find that four of the nine models pass the test (Table 2) with CanESM5 (Fig. A3), ACCESS-ESM1-5 (Fig. A4), and MIROC-ES2L (Fig. A5), showing a closer affinity to observations than the MPI-ESM1-2-LR model which we have focused on thus far.

A drawback of applying the KS test method here is that some model ensembles may be too small, and therefore, simply by chance, none of the ensemble members passes the KS test, even if the model is capable of producing such a simulation. Indeed, we see that four of the five models that are not passing the KS test have 10 or fewer historical ensemble members. For example, MPI-ESM1-2-HR is physically very similar to MPI-ESM1-2-LR, with the main difference being applying a higher resolution resulting in fewer simulated members, i.e., 10 instead of 30. Therefore, it is plausible that this model is also fit for purpose, but it was not run enough times to demonstrate that. The only exception to this rule is IPSL-CM6A, which, with 33 ensemble members, is more unlikely to have a plausible representation of the Amazonian vegetation dynamics.

Table 2 Observation percentiles p_cml in the mean τ_K model distribution, the fraction of ensemble members F_p>0.05 that scored p>0.05 in the KS test, and highest KS test p value of each ensemble. The higher the p value, the more likely the member resembles the observations. p_cml represents the fraction of the mean τ_K values of a model that are higher than the observational mean τ_K (see Fig. 4). Values are displayed for all nine historical ensembles and the control ensembles of the four well-agreeing historical ensembles. piControl ensembles are run under pre-industrial conditions. Ensembles that have good agreement have at least one member that does not significantly differ from the observations, according to the KS test.

Inspecting the Amazon mean trend of AR(1) for each of the ensemble members from the four models that passed the KS test, we find one model which displays a significant response to historical forcing (MIROC-ES2L; Fig. 4). To better show the shift between the pre-industrial control and historical simulations, we fit a bounded β distribution to the frequency distribution, and we use the α = 0.95 confidence interval to screen for significance. In all four models and for both simulations, the observed trend is on the edge of what is possible. In CanESM5, MPI-ESM1-2-LR, and ACCESS-ESM1-5, there is not a substantial difference between the forced and unforced simulations, but MIROC-ES2L shows a marked shift towards higher trends in the Amazon basin mean AR(1) when the model is exposed to historical boundary conditions.

We finally test the idea that an increasing AR(1) can be used as an early-warning signal of an eminent abrupt transition in the Amazon rainforest. This was done by inspecting the continuation of future projections of the two historical runs with largest positive and negative historical trends from MPI-ESM1-2-LR (Fig. 5). This is Shared Socioeconomic Pathway 5 with a radiative forcing of about 8.5 W m⁻² at the end of the century (SSP5-8.5). Under this strong forcing future scenario, the model exhibits a slightly increasing trend in AR(1), along with a decreasing trend of LAI of on average −0.05 m² m⁻². Neither of the two extreme runs nor any other ensemble members exhibit an abrupt decline in LAI. In fact, we do not observe a regime shift in any of the other four models that passed the KS test. The linear deterioration of forest viability can be caused by a variety of factors including limited water availability or progressing land use change (Fig. A7). However, the changes induced by future forcing do not facilitate the bifurcation-like behavior of the system.

4 Conclusions

In this study, we have tested the idea that trends in the persistence of vegetation density anomalies can be used as an early-warning signal for the Amazonian rainforest. This is particularly concerning against the backdrop of a large observed trend since 1991, suggesting that the forest has undergone a pronounced loss of resilience (Boulton et al., 2022). The trend in the anomaly persistence is quantified through the lag-1-year correlation, AR(1).

To this end, we inspect simulations from nine Earth system model ensembles initialized with different initial conditions in 1850, such that variations within each ensemble are an expression of the internal variability simulated by that particular model. We find that four of the models have ensemble members that are statistically indistinguishable from the observed trend. Of the other models with 10 or fewer realizations, 4 did not have a matching realization, and 1 model with 33 realizations clearly underperformed. Of the four well-performing model ensembles, three of them also showed trends similar to the observations in their unforced control simulations. These results suggest that the observed trend could simply be an expression of the internal variability and that longer data records would be needed to show that the opposite is the case.

This result is further corroborated by the spatial distribution of the increasing trend in AR(1) in the model simulations. Here it is found that ensemble members with substantial positive or negative trends show these in relatively large regions but not necessarily in those regions with large anthropogenic deforestation. This suggests that such anomalies could be associated with large-scale weather events.

We finally check whether trends in AR(1) can be used as an early-warning signal by investigating the relationship between the recent past and the rainforest evolution in future strong warming scenario projections (SSP5-8.5). However, there is no such relationship, and furthermore, none of the future simulations exhibits rapid transitions.

It is worth noting that even if the results presented here suggest that the Amazon rainforest has not lost its resilience and is unlikely to undergo bifurcations in the future, it does not mean the forest is invulnerable to human-caused stresses from global warming, deforestation, and fires. On the contrary, the results presented here suggest that global warming, as a major historical and future forcing in models, plays a minor role in the observed rapid decline in rainforest extent. Mitigation strategies to limit future rainforest loss could therefore be most effective when focusing on other, human-induced stress factors.

Appendix A

Figure A1 Kendall's τ_K distributions for three different VOD resolutions (0.25°×0.25° original VODCA). The means of these histograms are, respectively, ${\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{\mathrm{K},\mathrm{0.25}\mathrm{°}}$= 0.22, ${\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{K,\mathrm{0.5}\mathit{°}}$= 0.22, and ${\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{\mathrm{K},\mathrm{1}\mathrm{°}}$= 0.21. We use the 1°×1° resolution for our analysis since it is the closest to the model resolutions.

Figure A2 Maps of τ_K of MPI-ESM1-2-LR members, according to Fig. 3. The spatial mean for each run is indicated, and the observational value is ${\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{\mathrm{K}}$= 0.21 (resolution 1°×1°).

Figure A3 Cumulative distribution functions of τ_K for observations and CanESM5. Significant resemblance (p>0.05) according to the KS test is marked in red.

Figure A4 Cumulative distribution functions of τ_K for observations and ACCESS-ESM1-5. Significant resemblance (p>0.05) according to the KS test is marked in red.

Figure A5 Cumulative distribution functions of τ_K for observations and MIROC-ES2L. Significant resemblance (p>0.05) according to the KS test is marked in red.

Figure A6 LUCC difference between 2014 and 1990, as prescribed by LUH2 v2h (Hurtt et al., 2020), in CMIP6 models. (a) C3 annual crops (most small-seeded cereal crops), (b) primary forested land, (c) rangeland, and (d) pasture. Agricultural activity (a, c, d) increased particularly in the southeastern regions of the Amazon, while primary forested land (b) declined in the same areas.

Figure A7 LUCC difference between 2100 and 2015, as prescribed by LUH2 v2f SSP5 (Hurtt et al., 2020), in CMIP6 models. (a) C3 annual crops (most small-seeded cereal crops), (b) primary forested land, (c) rangeland, and (d) pasture. Particularly primary forested land reduction (b) effects increased in parts of the Amazon basin.

Figure A8 Anomaly of the spatially averaged AR(1) series for observations and the best-agreeing member for each historical ensemble. The anomaly is computed by subtracting the temporal mean from the time series and is plotted at the end of the 5-year sliding window. The best-agreeing member is shown in parentheses behind the model (the member number shown here agrees with the realization number of the CMIP6 variant ID) and chosen according to its trend statistic τ_K.