Drivers of and uncertainty in Amazon carbon sink long-term and interannual variability in CMIP6 models

The Amazon basin rainforest is a critical component of the climate system, currently representing 25 % of terrestrial carbon gains and storing 150×10⁹ to 200×10⁹ t of carbon. Whether the Amazon rainforest will remain a net carbon sink is an open scientific question: while its future stability and functioning may be compromised by climate change and anthropogenic pressures, Earth system model (ESM) divergence in the projections undermines the models' reliability in simulating its future evolution. In this study, we examined the contribution of different climatic drivers behind the long-term and interannual variability evolution of the carbon sink within the Amazon basin using 11 CMIP6 ESMs, shedding light on the main factors contributing to inter-model diversity. By adopting the carbon-cycle feedback framework with C4MIP experiments, our results underscore the dominant role of CO₂ fertilization in driving the long-term Amazon carbon sink trend and uncertainty. We also address the variability in carbon fluxes at the interannual timescale using a multivariate predictive model with historical and SSP5-8.5 ScenarioMIP simulations. In this respect, we emphasize the contribution of gross primary productivity (GPP) modulation by shortwave incoming radiation, which dominates net biome productivity (NBP) divergence across the ESM ensemble. Additionally, we demonstrate that temperature-driven anomalies will be the main mechanism responsible for the higher Amazon carbon sink sensitivity to the El Niño–Southern Oscillation (ENSO) under sustained global warming, predominantly as a result of the amplification of NBP sensitivity to temperature anomalies. Being that the representation of terrestrial carbon-cycle processes is still one of the main uncertainties undermining ESM projections, we therefore advocate for increased focus from modelling groups on a more accurate and consistent representation of land processes and parameterizations, which will hopefully lead to reduced uncertainties in simulations from the next generation of ESMs.

1 Introduction

The Amazon basin rainforest plays a fundamental role in the climate system, serving as a prominent actor in the land carbon cycle and exerting a significative influence on the global energy budget and hydrological cycle (Davidson et al., 2012). At the same time, the land carbon sink represents one of the crucial uncertainties affecting climate change future evolution (Canadell et al., 2021; Friedlingstein et al., 2006; Arora et al., 2020). Indeed, projections of the Amazon climate and carbon sink are still poorly constrained by state-of-the-art Earth system models (ESMs), indicative of persisting gaps in knowledge regarding this critical aspect of the Earth system (Ahlström et al., 2012; Baker et al., 2021; Koch et al., 2021a; Lin et al., 2023; Raoult et al., 2023; Xu et al., 2020; Zhu et al., 2019).

The Amazon ecosystem has been a long-term carbon sink during the past decades, contributing to approximately 25 % of terrestrial carbon gains, estimated at 0.42–0.65 GtC yr⁻¹ for the period of 1990–2007, a trend mainly driven since the 1980s by the CO₂ fertilization effect from rising CO₂ concentrations in the atmosphere (O'Sullivan et al., 2019; Pan et al., 2011; Phillips et al., 2009; Walker et al., 2021). Nevertheless, recent estimates have demonstrated a slowdown of net carbon sequestration and consequently a saturation and declining trend in the Amazon carbon sink, with an increase in carbon losses due to drought events and increased temperatures (Brienen et al., 2015; Hubau et al., 2020). Both CO₂ concentrations in the atmosphere and climatic conditions affect land carbon fluxes. Higher atmospheric CO₂ concentrations mainly exert a positive direct effect on photosynthesis through plants stomatal closure and the associated negative carbon–concentration feedback (Boer and Arora, 2009; Arora et al., 2020), while they can indirectly increase autotrophic respiration (R_a) and heterotrophic respiration (R_h; Gao et al., 2020). Climatic conditions, on the other hand, mainly affect vegetation carbon fluxes through temperature and water availability, with, for example, negative temperature impacts on primary productivity and a positive (negative) relationship between temperature (water availability) and changes in respiration rates (Canadell et al., 2021; Gentine et al., 2019; Green et al., 2019; Humphrey et al., 2018; Liu et al., 2020). In this respect, interannual variations in water availability and temperature in the Amazon basin are mainly related to the El Niño–Southern Oscillation (ENSO), which is responsible for a vast part of the climatological and net land CO₂ flux variability observed in tropical biomes (Bastos et al., 2018; Jones et al., 2001; Kim et al., 2016; Betts et al., 2020; Piao et al., 2020; Zhu et al., 2017). Accordingly, some of the most severe droughts observed in the Amazon basin in recent decades and the associated reduction in the net land CO₂ sink were forced by strong ENSO events (or El Niños), most prominently the 1997/1998 and 2015/2016 ones (Jiménez-Muñoz et al., 2016; Koren et al., 2018; Liu et al., 2017). Despite indications that temperature-driven gross primary productivity (GPP) anomalies were responsible for a decreased Amazonian carbon sink in the 2015/2016 event (Bastos et al., 2018; Zhang et al., 2019), it is still currently debated whether fluctuations in temperatures or water availability are the dominant drivers of interannual carbon variability in tropical biomes, with recent research indicating the increased importance of water availability as a controlling factor in the past decades and with ESMs failing to reproduce this observed behaviour (Humphrey et al., 2018; Jung et al., 2017; Liu et al., 2023; Zhang et al., 2023).

Consequently, at least three factors will contribute to determining whether, and to what extent, the Amazon ecosystem will remain a net carbon sink in future decades under sustained positive radiative forcing: mean-state climatic changes, nutrient limitation, and ENSO modulation. First, a significant increase in surface air temperature and a marked decline in water availability in the Amazon basin deriving from scenarios of increased greenhouse gas emission will most likely result in a less effective carbon sink by the end of the 21st century. In particular, coupled climate models suggest that the reduction in precipitation is driven by reduced evapotranspiration resulting from stomatal closure response to increased CO₂, which leads to changes in local surface energy balance and atmospheric circulation patterns (Kimm et al., 2024; Kooperman et al., 2018; Langenbrunner et al., 2019; Li et al., 2023). Then, nitrogen and phosphorous will likely limit tropical forest productivity (Fleischer et al., 2019), an effect partly counterbalanced by the positive influence exerted by the increasing atmospheric CO₂ concentrations (Huntingford et al., 2013; Koch et al., 2021b). Lastly, an increased Amazon vegetation sensitivity to ENSO is expected under a range of global warming scenarios (Kim et al., 2017; Park et al., 2020; Uribe et al., 2023). In particular, global warming could affect the ENSO–Amazon teleconnection either from changes in the mean state and extremes in the ENSO cycle or from variations in the ENSO teleconnection mechanism itself (Beobide-Arsuaga et al., 2021; Cai et al., 2021; Chen et al., 2017; Taschetto et al., 2020; Yeh et al., 2018; Zheng et al., 2017). Notably, ENSO amplitude is slightly yet significantly enhanced under future global warming scenarios (Beobide-Arsuaga et al., 2021; Cai et al., 2022), and regional patterns of precipitation and temperature anomalies over South America associated with ENSO teleconnections are projected to be amplified in warmer climates (McGregor et al., 2022; Parsons, 2020; Perry et al., 2020;).

Given these premises, in this research we investigate the Amazon carbon sink using coupled simulations from CMIP6 ESM simulations (Eyring et al., 2016; O'Neill et al., 2016) to separate the relative contributions of long-term changes and interannual variability under sustained atmospheric CO₂ concentrations and global warming. Specifically, throughout the paper, two questions that remain underexplored in the literature are tackled: what are the relative contributions of interannual variability (IAV) and long-term climatic changes in the projected Amazon net carbon sink evolution? What is the role of ENSO in exerting a control on temperature and water availability? In tackling these questions, an attempt is made to identify the factors contributing to inter-model diversity in Amazon vegetation productivity projections across the ESMs considered.

2 Data and methods

2.1 Data

We use simulations from 11 CMIP6-generation ESMs that contributed to both C4MIP and ScenarioMIP projects (Eyring et al., 2016; Jones et al., 2016; O'Neill et al., 2016). In particular, the following experiments have been considered: 1pctCO2-bgc, 1pctCO2-rad, ssp585-bgc, ssp585-rad, historical and SSP5-8.5. ScenarioMIP simulations (historical and SSP5-8.5) aim to reproduce the climatic response to realistic forcing of the historical period and to a prescribed 8.5 W m⁻² radiative forcing increase by the end of the 21st century. C4MIP experiments (1pctCO2-bgc, 1pctCO2-rad, ssp585-bgc, ssp585-rad) are idealized concentration-driven carbon–climate simulations generated to better understand and quantify changes in the ocean and land carbon storage and fluxes under different climatic conditions. Specifically, the experiments aim to test the carbon-cycle response to the effect of increased CO₂ concentrations in the atmosphere (simulations with the suffix “-bgc”) and increased radiative forcing resulting from higher atmospheric CO₂ concentrations affecting the climate system (simulations with the suffix “-rad”). These two categories of simulations differ in their model set-up so that, in the former (-bgc), only the model land and ocean carbon cycle respond to a CO₂ increase, while the radiation scheme uses a preindustrial CO₂ concentration, therefore allowing the testing of the biogeochemical effect of atmospheric carbon dioxide increase without the associated radiative forcing. Reversely, in -rad simulations the biogeochemical effect is factored out; hence the climate responds to the radiative forcing by increased CO₂ concentration, whereas the carbon cycle remains constrained by a preindustrial atmospheric CO₂ level. In our analysis, we adopt C4MIP simulations forced either with a 1 % yr⁻¹ increase in atmospheric CO₂ concentrations up to 4 times the preindustrial level of 280 ppm, with no confounding effect of changes in land use, non-CO₂ greenhouse gases and aerosols, or with a standard CO₂ pathway from the SSP5-8.5 scenario. For the sake of our goal, we consider these differences negligible, as considering both 1pctCO₂ and SSP5-8.5 experiments allows us to have a higher ensemble of data available.

Those ESMs for which at least one realization is available for all the simulation experiments and that have all the prognostic variables needed for the analysis available (or computable indirectly) have been included in the research. Table 1 presents an overview of the experiments considered and the methodologies adopted for each experiment. For the historical and SSP5-8.5 experiments, specifically, the first five realizations of each model (when more than one was available) were used to account for the uncertainty stemming from internal climatic variability, with the limit of having the same simulation members for the historical and SSP5-8.5 scenarios to make a pairwise comparison. The details of the ESMs used are reported in Table S1. For these historical and SSP5-8.5 simulations, the analyses have been performed on single realizations and aggregated by model solely for the graphical presentation of the results. The evolution of the land carbon cycle is investigated by considering long-term changes and interannual variability in net biome productivity (NBP) values, which represent the balance of gross primary productivity (GPP) due to photosynthesis at the net of autotrophic respiration (R_a), heterotrophic respiration (R_h) and disturbances, such as fire dynamics (for those model having a fire module) and land use change (LUC). Overall, the following variables have been considered in the study: sea-surface temperature (SST), net biome productivity (NBP), gross primary production (GPP), autotrophic respiration (R_a), heterotrophic respiration (R_h), precipitation (Pr), soil moisture (mrso), air surface temperature (T) and shortwave incoming radiation (SW_in).

Table 1 Overview of the experiments and methodologies used in the research. We refer to the work of Jones et al. (2016) and O'Neill et al. (2016) for more information on the C4MIP and ScenarioMIP experiments reported in this table.

Zonal statistics computed based on measurements made within the Amazon basin and presented throughout the results are obtained by considering the grid cells confined within the Amazon basin shapefile, available from the SO HYBAM service (INPE, 2019, https://hybam.obs-mip.fr/, last access: 17 March 2024).

The ESM performances in representing ENSO properties, the Amazon climatology, and carbon and energy fluxes are evaluated against observational and quasi-observational products and are reported in Figs. S1–4. The HadISST dataset is used for assessing ESM sea-surface temperatures (Rayner et al., 2003), while ERA5 and ERA5-Land products are used to validate temperature, precipitation and soil moisture (Hersbach et al., 2020; Muñoz-Sabater et al., 2021). Finally, the FLUXCOM-RS + METEO dataset, specifically the one forced with the WFDEI meteorological dataset (Weedon et al., 2014), is used as a reference for both carbon fluxes (GPP; net ecosystem productivity, NEP; total ecosystem respiration, TER) and energy fluxes (shortwave incoming radiation) (Jung et al., 2019, 2020).

To evaluate ESMs against the reanalysis products, as well as to compute multi-model means across variables, a conservative remapping algorithm is applied to all the data to get a regular 1° longitude × 1° latitude grid, with the exception of SST from ESMs with a curvilinear grid (Table S1), for which a distance-weighted (nearest-neighbour) average remapping is applied. The validation procedure refers to the climatological period of 1979–2013. When comparing the carbon fluxes from ESMs with FLUXCOM data, the total ecosystem respiration is obtained by summing the contributions of R_a and R_h. An overview of the ESM evaluation performances is available in the Supplement (Figs. S1–4).

2.2 Long-term mean-state climatic effects

Increasing trends of atmospheric CO₂ concentrations affect terrestrial carbon sinks directly through a fertilization effect on vegetation (carbon–concentration feedback) and indirectly by forcing changes in the physical climate via a strengthened greenhouse effect, which in turn affects vegetation (carbon–climate feedback). A common approach to disentangle the two effects relies on the carbon-cycle feedback framework, by which it is possible to estimate the magnitude of the carbon–concentration feedback and the carbon–climate feedback (Jones et al., 2016; Arora et al., 2020). The contribution of these two feedbacks is estimated with the following linear equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\beta }={\displaystyle \frac{\mathrm{\Delta }{\mathrm{NBP}}_{\mathrm{cum},\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{ppm}}_{\mathrm{BGC}}}},\text{(2)}& {\displaystyle }\begin{array}{rl}& {\mathit{\gamma }}_T={\displaystyle \frac{\mathrm{\Delta }{\mathrm{NBP}}_{\mathrm{cum},\mathrm{RAD}}}{\mathrm{\Delta }{T}_{\mathrm{RAD}}}},{\mathit{\gamma }}_{\mathrm{mrso}}={\displaystyle \frac{\mathrm{\Delta }{\mathrm{NBP}}_{\mathrm{cum},\mathrm{RAD}}}{\mathrm{\Delta }{\mathrm{mrso}}_{\mathrm{RAD}}}},\\ & {\mathit{\gamma }}_{{\mathrm{SW}}_{\mathrm{in}}}={\displaystyle \frac{\mathrm{\Delta }{\mathrm{NBP}}_{\mathrm{cum},\mathrm{RAD}}}{\mathrm{\Delta }{\mathrm{SWin}}_{\mathrm{RAD}}}}.\end{array}\end{array}$$

The carbon sink is represented here by cumulative NBP, whose long-term sensitivity to CO₂ (ppmv; β), and climate (γ) is estimated from biogeochemical-only coupled simulations (1pctCO₂-bgc and ssp585-bgc) and radiative-only coupled simulations (1pctCO₂-rad and ssp585-rad), respectively. Within the former, the increased CO₂ atmospheric concentration uniquely exerts a biogeochemical effect which impacts the terrestrial and ocean carbon-cycle processes, whereas in the radiative simulation, only the radiative transfer processes in the atmosphere are affected by the changing atmospheric CO₂ concentrations, with no consequences for biochemical processes (Friedlingstein et al., 2006; Jones et al., 2016). Changes are expressed with respect to the first year of the simulations; for 1pctCO₂ experiments, only the years with atmospheric CO₂ ranges similar to the SSP5-8.5 scenario (400–1135 ppmv, resulting in 104 years) have been considered. Lastly, the carbon–climate feedback has been computed with respect to surface atmospheric temperatures (T), whereas the carbon sink sensitivity to soil moisture (mrso) and shortwave incoming radiation (SW_in) is used to represent the variety of mean-state climatic changes affecting the cumulative carbon fluxes within the Amazon basin.

2.3 Carbon fluxes sensitivity at interannual timescales

We further assessed the sensitivity of carbon fluxes within the Amazon basin at interannual timescales to explore the short-term variability contribution of different climatic factors. In this respect, two additional analyses have been performed. A first one aimed at assessing the relative contribution of temperature, soil moisture and shortwave incoming radiation on NBP, GPP, R_a and R_h by adopting a multivariate predictive regression model as in the following equation, for every ESM considered:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Delta }\mathrm{NBP}& ={\displaystyle \frac{\mathit{\delta }\mathrm{NBP}}{\mathit{\delta }T}}\cdot \mathrm{\Delta }T+{\displaystyle \frac{\mathit{\delta }\mathrm{NBP}}{\mathit{\delta }\mathrm{mrso}}}\cdot \mathrm{\Delta }\mathrm{mrso}\\ & +{\displaystyle \frac{\mathit{\delta }\mathrm{NBP}}{\mathit{\delta }{\mathrm{SW}}_{\mathrm{in}}}}\cdot \mathrm{\Delta }{\mathrm{SW}}_{\mathrm{in}}+\mathit{ϵ}.\end{array}\end{array}$$

Temperature, soil moisture and shortwave incoming radiation have been chosen because of their primary control on carbon flux interannual variability (Jung et al., 2017). We used soil moisture rather than precipitation as a proxy for water availability due to its stronger control on terrestrial carbon fluxes (Humphrey et al., 2018). We acknowledge that latent heat fluxes could partially contribute to interannual NBP variability by influencing the vapour pressure deficit (VPD), which in turn affects NBP. However, we decided to exclude this factor from our framework (Eq. 3) because transpiration, which is regulated by stomatal conductance, impacts both NBP and latent heat fluxes, likely introducing correlation between them.

All the variables considered have been averaged over the calendar year, as the focus of the analysis is on interannual variability.

The dependent variable (NBP, as well as GPP, R_a and R_h) has been linearly detrended to remove the influence of CO₂ fertilization and other long-term climatic changes, while the three independent variables have been detrended and standardized (subtracting the mean and dividing by their standard deviation) to obtain comparable coefficients. Therefore, the standardized coefficients represent the quantitative contribution of temperature, soil moisture and shortwave incoming radiation considering the simultaneous confounding impacts of the other variables. Specifically, a 5-fold cross-validation ridge regression model has been adopted, due to the penalty score in the cost function of ridge regression that helps to account for the collinearity among the predictors themselves. The cross-validation procedure entailed a randomized split of the dataset into a training set (80 %) and validation set (20 %) and allowed for the selection of the best-performing regularization parameter among a set of values spaced evenly on a log scale. The multivariate model is derived from the scikit-learn package, available for the Python scripting language (Pedregosa et al., 2012).

Given the predominant modulation of carbon flux IAV within the Amazon basin by means of ENSO (Betts et al., 2020), a further analysis is conducted to assess the ENSO contribution to Amazon vegetation productivity mediated by either temperature or water availability (using soil moisture as a proxy). For this, an annual time series of ENSO is obtained for each historical and SSP5-8.5 realization by averaging the corresponding monthly Nino3.4 index over the calendar year. The Nino3.4 index is defined as the 5-month moving average of mean sea-surface temperatures over the region 170–120° W and 5° S–5° N, subsequently detrended by means of a first-order polynomial to remove the warming trend of SST. By taking the univariate sensitivities of temperature and soil moisture to ENSO and accounting for the partial derivatives of NBP with respect to temperature and soil moisture as in Eq. (3), we estimate the contribution of ENSO driven by the two mechanisms as below:

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\delta }{\mathrm{NBP}}_{\mathrm{n}\mathrm{34}}^T={\displaystyle \frac{\mathit{\delta }\mathrm{NBP}}{\mathit{\delta }T}}\cdot {\displaystyle \frac{\mathrm{d}T}{\mathrm{dn}\mathrm{34}}},\text{(5)}& {\displaystyle }\mathit{\delta }{\mathrm{NBP}}_{\mathrm{n}\mathrm{34}}^{\mathrm{mrso}}={\displaystyle \frac{\mathit{\delta }\mathrm{NBP}}{\mathit{\delta }\mathrm{mrso}}}\cdot {\displaystyle \frac{\mathrm{dmrso}}{\mathrm{dn}\mathrm{34}}}.\end{array}$$

Using the partial derivatives $\frac{\mathit{\delta }\mathrm{NBP}}{\mathit{\delta }T}$ and $\frac{\mathit{\delta }\mathrm{NBP}}{\mathit{\delta }\mathrm{mrso}}$ ensures that the effects of temperature and soil moisture are considered independently, accounting for the potential confounding influences of one on the other, a condition that would not be met if the univariate estimates $\frac{\mathrm{dNBP}}{\mathrm{d}T}$ and $\frac{\mathrm{dNBP}}{\mathrm{dmrso}}$ were applied. A Mann–Whitney U test of independence with Bonferroni correction was used to assess whether the zonal values of the regression coefficients within the Amazon basin are significantly different between the historical period and SSP5-8.5 scenario. To mitigate the risk of overstating the significance of the statistical tests conducted, we employ a false discovery rate (FDR) control method based on Wilks (2016). This approach effectively addresses the issue of multiple hypothesis testing, ensuring a more accurate interpretation of the obtained results.

3 Results and discussion

3.1 Inter-model uncertainties in carbon and climatic drivers

Amazon basin vegetation productivity and climatology show substantial overlap across models during the historical period but strongly diverging trends across models during the SSP5-8.5 scenario, with individual models differing in the magnitude and sign of projected changes (Fig. 1a). The multi-model ensemble yields a cumulative NBP mean by the end of the 21st century of 59.3 ± 62 gC m⁻². Considering the projections of the cumulative carbon sink in Fig. 1a and carbon fluxes in Fig. S5, the inter-model uncertainty is much higher compared to intra-model uncertainty (which stems from the intrinsic climatic variability expressed in each realization and is represented by the ± 1 standard deviation spread in Figs. 1 and S5). For the physical variables in Fig. 1b–e, intra-model uncertainty is considerably higher than for carbon fluxes and reflects the substantial internal climate variability intrinsic to each simulation. These considerations already highlight that part of the divergence across carbon-cycle predictions is related to differences in the land sensitivity to climatic forcings, rather than to uncertainties in the evolution of the climate itself. Overall, the climatological variables present stronger agreement and coherence in the sign of projected changes among ESMs with respect to NBP. A multi-model mean reduction of −21.13 mm per month and −128.93 kg m⁻² is projected for precipitation and soil moisture, respectively (Fig. 1c, d), while an increase is observed for temperature and incoming shortwave radiation (+7.08 °C and +9.23 W m⁻², respectively; Fig. 1b, e). Still, the multi-model ensemble spread at the end of the 21st century remains substantial for all these variables. Differences among models of 1 or even 2 orders of magnitude could be seen for instance between MRI-ES2L and CanESM5 for shortwave incoming radiation or for MIROC-ES2L and IPSL-CM6A-LR regarding soil-moisture projections.

For some models, divergent NBP projections cannot be easily attributed to similar deviations in the individual climatic drivers. For instance, ACCESS-ESM1-5 and NorESM2-LM project similar end-of-the-century temperature and soil-moisture changes but strongly diverging trends of cumulative NBP.

Figure 1 Long-term trends of (a) cumulative NBP, (b) temperature, (c) precipitation, (d) soil moisture and (e) shortwave incoming radiation in the Amazon basin for the historical and SSP5-8.5 experiments. Trends are computed with respect to the first 30-year mean of the historical period (1850–1880) and are visualized as a 10-year moving average for clarity. For the models with more than one realization, both the model ensemble mean (line) and the spread (± 1 standard deviation, shading) are shown.

Regarding the intra-model spread, the highest influence of internal climatic variability (± 1 standard deviation) is observed in precipitation and shortwave incoming radiation, followed by soil moisture (spreads in Fig. 1c, e and d). This indicates that within the Amazon basin, the major source of uncertainty deriving from internal climate variability is associated with cloud formation and coverage, which is causally linked with the amount of precipitation (thus soil-moisture content) and shortwave incoming radiation within the regional domain. All the carbon fluxes from which NBP is derived (GPP, R_a and R_h) depict an increasing trend throughout the 21st century (Fig. S5), with the notable exception of GPP and heterotrophic respiration for ACCESS-ESM1-5. Among these carbon fluxes, R_h presents the highest end-of-21st-century normalized uncertainty (186.74 gC m⁻², z score SD of 1.62), followed by GPP and R_a (753.59 and 548.71 gC m⁻², z score SD of 1.12 and 0.76, respectively). This shows that uncertainty in cumulative NBP does not solely stem from uncertainty in single climatic factors; instead, it arises from inconsistencies and differences in how models represent photosynthetic activity and autotrophic and heterotrophic respiration. Overall, inconsistencies in projected Amazon NBP cannot be simply understood as a consequence of discrepancies in the climatic factors affecting vegetation, as discussed already by Heavens et al. (2013). A point that deserves attention is therefore how the projected carbon sink in the Amazon basin is sensitive to mean-state changes in environmental drivers and climatic variability at the interannual timescale.

3.2 Long-term carbon sink sensitivity

By applying the carbon-cycle feedback framework (see “Data and methods”), the net carbon sink (cumulative NBP) response is explained by the factor β (CO₂ fertilization) and by the factor γ (climate effect). The net carbon sink trend driven uniquely by the CO₂ fertilization effect and climate is reported in Fig. S6a and b, respectively, where the consistent agreement across ESMs on the direction of the projected land carbon response to these two factors is evident. The spatial magnitude of β derived from the 1pctCO₂-bgc simulation is shown in Fig. S7. Most of the ESMs project a strong CO₂ fertilization effect within the Amazon basin for the range of atmospheric CO₂ concentrations of the SSP5-8.5 scenario, despite considerable variability that could be observed across the models. CESM2 and NorESM2-LM present the higher vegetation sensitivity to atmospheric CO₂, which is reasonable due to the fact that both models share the same land module, CLM5 (Table S1). CMCC-ESM2 also presents high β values, but these are restricted to the most forest-dense pixel cells, as is clear from the sharp decline in the southern and south-eastern part of the basin, a pattern which is even more drastic in ACCESS-ESM1-5.

The carbon-cycle feedback framework entails describing the positive carbon–climate feedback uniquely considering surface air temperature (γ_T), which is therefore considered an overall representation of long-term climate impacts, as shown in Fig. 2. We additionally focus on different explainable variables as equivalent terms representing long-term climate impacts (γ_mrso and γ_SWin). Thus, despite γ_T, γ_mrso and γ_SWin not being intended as cumulative long-term impact coefficients, the values of their standardized coefficients are reported in Fig. S8 with the purpose of providing a quantitative comparison of the ESM sensitivity to diverse climatic factors. Remarkable negative effects on the Amazon carbon sink are associated with long-term increasing temperatures and partially, with less intensity, with rising trends in shortwave incoming radiation. The positive co-variability in the carbon sink with soil moisture is almost the mirror of the cumulative NBP sensitivity to shortwave incoming radiation, a consequence of the strong inverse correlation among the two factors and reflecting a reduction in cloud coverage and precipitation in the Amazon basin with increasing atmospheric CO₂ and climate warming.

Consequently, the long-term cumulative NBP for the 1pctCO₂ scenario is reconstructed by means of its sensitivity to mean-state changes in temperature and atmospheric CO₂ and is reported in Fig. 2. These results emphasize that CO₂ fertilization effects on vegetation productivity are expected to dominate and overcome the negative temperature influence on the long-term carbon sink evolution in the SSP5-8.5 scenario.

Figure 2 Contribution of carbon–concentration (β) and carbon–climate (γ_T) feedbacks to net carbon sink projections at the end of the 1pctCO2-bgc and 1pctCO2-rad simulations. Black dots represent the net cumulative NBP resulting from the arithmetic sum of cumulative NBP from the 1pctCO₂-bgc and 1pctCO₂-rad simulations, whereas the blue crosses are the result of the arithmetic sum of the estimated cumulative NBP by adopting the univariate β and γ_T coefficients. Values are averaged within the Amazon basin domain.

Regarding the differences across ESMs, these discrepancies in the long-term climate-driven trend of the net carbon sink in the Amazon region emerge from a very strong correlation between its sensitivity to soil moisture γ_mrso and temperature γ_T, as reported in Fig. 3 and Fig. S9 for both the 1pctCO2-rad and the ssp585-rad simulations. Therefore, the Amazon carbon sensitivity to long-term climatic changes not only is discernible in the negative temperature effect, but also is found in the positive influence on soil moisture, indicating that the land modules of the ESMs are similarly sensitive to both the factors, despite having an inverse relationship. Additionally, temperature emerges as playing a slightly stronger role in describing the variations in the carbon sink response across the different ESMs compared to soil moisture (R² coefficient of 0.85 vs. 0.71, as illustrated in Fig. 3a and b).

Figure 3 Inter-model uncertainty in climate-driven cumulative NBP (y axis) as explained by differences in ESM representation of temperature (γ_tas) and soil moisture (γ_mrso) impacts in 1pctCO2-rad and ssp585-rad simulations. The reported values are zonally averaged within the Amazon basin geographical domain.

3.3 Interannual variability in carbon fluxes

Amazon carbon sink variability on shorter (interannual) timescales is assessed for the historical and SSP5-8.5 scenarios considering the modulation of temperature, soil moisture and shortwave incoming radiation anomalies using a multilinear ridge regression predictive model (see “Data and methods”). The skills of the model optimized with the 5-fold cross-validation procedure are reported in Fig. S10. Overall, a multi-model mean coefficient of determination of 0.55 is obtained, despite substantial differences across ESMs. Remarkably lower values of explained variance (R²<0.4) are found for CESM2, CMCC-ESM2 and NorESM2-LM, whereas CanESM5, MIROC-ES2L and UKESM1-0-LL stand out as the models with the highest goodness of fit and predictive capability (R²>0.6) (Fig. S10a). Undoubtedly, we acknowledge the limitations that arise from the adoption of a linear ridge regression model. As we did not account for interactive and non-linear effects among the predictors influencing NBP IAV, we are not able to capture a portion of the unexplained variance in our regression model. The results in Fig. S10 reflect this fact, suggesting that other factors, as well as the effects emerging from the interaction of temperature, soil moisture and shortwave incoming radiation, likely have an important influence on NBP IAV, especially for CESM2, NorESM2-LM and partly CMCC-ESM2. Nevertheless, we opted for this modelling framework, as it still allows our results to be compared with prior works (Humphrey et al., 2018; Jung et al., 2017).

Figure 4 Partial derivatives explaining the contribution of temperature (a), soil moisture (b) and shortwave incoming radiation (c) to interannual NBP, averaged across the Amazon basin. The black vertical bars represent the spread in the predictor coefficients for models with more than one realization available, whereas the asterisks indicate the level of significance (p value), averaged over the Amazon basin, associated with every coefficient. Statistical significance refers to the following convention: ^* $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{2}}<p\le \mathrm{5.00}\times {\mathrm{10}}^{-\mathrm{2}}$; ${}^{**}$ $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{3}}<p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{2}}$; ${}^{***}$ $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{4}}<p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{3}}$; ${}^{****}$ $p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{4}}$.

The relative importance of the three variables in explaining the NBP interannual variability is reported in Fig. S11 for both the historical period (panel a) and the SSP5-8.5 scenario (panel b), revealing the dominant regulation of temperature for all the ESMs except CanESM5, IPSL-CM6A-LR and MPI-ESM1-2-LR, which primarily have a soil-moisture-driven variability, and MIROC-ES2L, which stands out for its particularly high variance explained by shortwave incoming radiation and its low contribution of soil moisture. As will be discussed later in the next section, temperature modulation is expected to be particularly predominant in the SSP5-8.5 scenario (Fig. S11b), intensified by climate change. A similar conclusion can be drawn by observing the standardized coefficients estimated for NBP. Those, reported for the SSP5-8.5 scenario in Fig. 4, point to inter-model qualitative agreement with respect to the contribution of temperature (negative effect) and soil moisture (positive effect), despite remarkable differences in the magnitude of the drivers' contributions. Lower net carbon sink variability, on the other hand, is associated with shortwave incoming radiation, with ESMs diverging on the sign of the partial derivative (Fig. 4c).

Similar results are found for the multilinear regression estimation in the historical period; compared to this, a general increase in the sensitivity is found for temperature and soil-moisture effects in the SSP5-8.5 scenario (Fig. S12), which is mainly a consequence of increased future interannual NBP variability. In Fig. 4, it can be seen that 7 out of 11 models project temperature as the first NBP predictor, while for CanESM5, IPSL-CM6A-LR and MPI-ESM1-2-LR, soil-moisture results are more important and MIROC-ES2L presents a very similar contribution from both temperature and shortwave incoming radiation. Overall, these results complement and partially confirm what has been shown in recent research (Padrón et al., 2022), which found temperature to be more important than soil moisture in explaining interannual variability in NBP in an ensemble of models, considering the SSP1-2.6 low-emission scenario.

The spatially explicit multi-model ensemble coefficients for NBP are reported in Fig. 5, with hatching delimiting those pixel cells where 8 out of 11 ESMs agree in the sign of the partial derivative. Temperature impacts are particularly evident in the northern part of the Amazon, whereas the positive influence of soil moisture presents high agreement almost everywhere, with the exception of the north-western part of the basin. Considering the almost specular multi-model agreement in the positive value of shortwave incoming radiation coefficients in that specific domain, this could reflect the fact that the majority of the ESMs present an energy-limited regime in the area, with shortwave incoming radiation (rather than water availability) being the main limiting factor constraining ecosystem productivity at interannual timescales.

Figure 5 Multi-model ensemble mean of the coefficient values for the climatic drivers obtained by the multilinear regression, for the SSP5-8.5 scenario. Hatching represents those grid cells for which at least 8 out of 11 ESMs agree in the sign of the predictor value. The Amazon basin, obtained from the SO HYBAM service (https://hybam.obs-mip.fr/, last access: 17 March 2024), is also represented.

As observed for long-term effects in Fig. 3, the inter-model differences in the projection of net carbon sink variability at interannual timescales are characterized by an inverse relationship between temperature and soil-moisture influence: high NBP interannual variability (expressed as NBP standard deviation, the y axis in the panels of Fig. S13) is driven by lower temperature (Fig. S13a) and higher soil-moisture contributions (Fig. S13b).

We then identify the physical processes associated with the net carbon sink that exhibits the greatest uncertainty by analysing the variability in the estimated coefficients for each dependent variable (NBP, GPP, R_a and R_h), reported in Table 2.

In this procedure, the coefficients obtained from the multivariate regression model have been further standardized by dividing each coefficient by the standard deviation of its corresponding dependent variable, ensuring a fair comparison both across models and for the different dependent variables.

Table 2 Inter-model uncertainty within the Amazon basin, expressed as the standard deviation of the partial derivative values across all the ESMs. Rows represent the partial derivatives of temperature, soil moisture and shortwave incoming radiation, with respect to the dependent variables reported in the columns (carbon fluxes).

The high disagreement in the sign of the shortwave incoming radiation coefficients found for NBP (panel c in Fig. 4) originates from GPP sensitivity to shortwave incoming radiation $\frac{\mathit{\delta }\mathrm{GPP}}{\mathit{\delta }{\mathrm{SW}}_{\mathrm{in}}}$ (Fig. S14). On the one hand, uncertainties in net carbon sink projections across ESMs therefore arise, mainly from differences in photosynthesis modulation by light availability. On the other hand, heterotrophic respiration sensitivity to soil moisture and temperature represents another prominent source of uncertainty across ESMs (see also panels a and b in Fig. S15 showing the disagreement across ESMs coefficients), whereas autotrophic respiration shows median contributions to net carbon sink uncertainty (Table 2 and Fig. S16).

Figure 6 Variability in NBP at interannual timescales associated with ENSO (y axis) and mediated by either soil moisture or surface air temperature (x axis), for every ESM. Reported is the distribution of the values within the Amazon basin, for the historical period (light green) and the SSP5-8.5 scenario (orange). Statistical differences among the distribution of the coefficients for the two periods are tested by means of a Mann–Whitney U test and are reported by the asterisks above the plots with the following convention: ^* $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{2}}<p\le \mathrm{5.00}\times {\mathrm{10}}^{-\mathrm{2}}$; ${}^{**}$ $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{3}}<p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{2}}$; ${}^{***}$ $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{4}}<p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{3}}$; ${}^{****}$ $p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{4}}$.

Figure 7 Multi-model ensemble mean of (a) temperature-mediated and (b) soil-moisture-mediated ENSO influence on NBP interannual variability. Only grid cells within the Amazon basin are included. Statistical significance, tested by means of a Mann–Whitney U test, is reported in the asterisks above the plot, and refers to the following convention: ^* $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{2}}<p\le \mathrm{5.00}\times {\mathrm{10}}^{-\mathrm{2}}$; ${}^{**}$ $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{3}}<p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{2}}$; ${}^{***}$ $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{4}}<p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{3}}$; ${}^{****}$ $p\le \mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{4}}$.

3.4 ENSO modulation

Given the predominant role of large-scale climatic modes of variability in influencing the local climate and carbon cycle, the NBP response to ENSO has been investigated. ENSO is the most important climatic mode affecting and modulating net tropical carbon sink variability at interannual timescales, and changes in key properties of ENSO could significantly impact the Amazon basin region. For example, an increase in the frequency or in the magnitude of El Niño and La Niña events under global warming, a possibility explored in previous studies (Berner et al., 2020; Brown et al., 2020; Cai et al., 2014, 2015; Fredriksen et al., 2020), could have important implications, especially for the Amazon ecosystem due to a stronger inhibition of the boreal winter southward shift of the Intertropical Convergence Zone (ITCZ) that regulates the ENSO tropical teleconnection pathway. Another ENSO property that could influence the Amazon carbon cycle is the spatial location of SST anomalies in the tropical Pacific, resulting in the distinction between Central Pacific (CP) and Eastern Pacific (EP) ENSO events (Capotondi et al., 2020). This spatial diversity in ENSO has been shown to have different impacts on GPP and NEP in many regions of the world, especially in the Amazon basin (Dannenberg et al., 2021), and some studies have also shown that CP El Niños are projected to be more frequent under 21st-century warming (e.g. Shin et al., 2022). However, the consideration of these ENSO properties lies outside the domain and the purpose of the present research, and we leave it for future studies.

The ESMs in the considered ensemble largely overestimate the observed ENSO amplitude in the interannual-to-decadal band, with the associated spectra typically featuring a narrow peak around the 3-year periodicity (Fig. S17). The ESMs also yield a diversity of changes in ENSO spectral characteristics under the warming scenario: generally, all the models represent a shift of the ENSO signal towards higher frequencies, with weaker amplitude at decadal timescales and stronger amplitude at interannual timescales (Fig. S17). Changes in the amplitude of the ENSO signal, represented by the Nino3.4 index standard deviation, between the SSP5-8.5 scenario (orange dots) and the historical period (light-green dots) are shown in Fig. S18; 9 out of 11 ESMs show an increased ENSO variability under the SSP5-8.5 scenario, whereas CESM2 does not project any relevant changes and UKESM1-0-LL predicts a slight decrease in ENSO amplitude.

We estimate the NBP variability related to ENSO by considering both the sensitivity of NBP to temperature and soil moisture and the sensitivity of temperature and soil moisture to the Nino3.4 signal (see “Data and methods”), thus allowing the decomposition of the impacts associated with the two different mechanisms. The results are shown in Fig. 6. First, all the models represent larger carbon sink losses related to temperature-driven anomalies compared to soil-moisture ones. Second, temperature and soil-moisture impacts driven by ENSO are expected to increase in magnitude during the future SSP5-8.5 scenario compared to the historical period, with the increase in $\mathit{\delta }{\mathrm{NBP}}_{\mathrm{n}\mathrm{34}}^T$ appearing to be statistically significant for all the models but MPI-ESM1-2-LR and NorESM2-LM. Regarding $\mathit{\delta }{\mathrm{NBP}}_{\mathrm{n}\mathrm{34}}^{\mathrm{mrso}}$, on the contrary, most of the ESMs (7 out of 11) project either a non-significant change or a relatively low statistical significance (0.05<p value <0.01). While the carbon-cycle response to ENSO within the Amazon region presented here is well in agreement with previous research (Betts et al., 2020; Kim et al., 2017; Le, 2023; Le et al., 2021), our results further identify temperature as the key factor in the mechanism by which ENSO affects Amazon carbon fluxes in a high-radiative-forcing future scenario. Specifically, the higher temperature-mediated ENSO impacts arise from the higher (lower) co-variability observed in the ESMs between temperature (soil-moisture) anomalies and ENSO, both during the historical period and in SSP5-8.5 projections (Fig. S19). This is accentuated by the stronger response of NBP to temperature changes compared to soil moisture for the majority of ESMs, as shown in Fig. 4.

The influence of ENSO by means of temperature and soil-moisture anomalies is additionally examined considering the multi-model ensemble mean of $\mathit{\delta }{\mathrm{NBP}}_{\mathrm{n}\mathrm{34}}^T$ and $\mathit{\delta }{\mathrm{NBP}}_{\mathrm{n}\mathrm{34}}^{\mathrm{mrso}}$. To do so, the ESMs coefficients reported in Fig. 4 have been re-gridded to a common 1×1 grid using bilinear interpolation. As discernible from Fig. 7, the temperature-driven impacts of ENSO are expected to be clearly stronger in the future SSP5-8.5 scenario compared to the historical period (panel a), as well as with respect to soil-moisture-driven impacts (panel b). Considering that climatic variability associated with ENSO in the Amazon basin is not expected to significantly change in the SSP5-8.5 scenario with respect to the historical simulation for all the models (Fig. S19), these results indicate that climate change will significantly amplify the NBP sensitivity to temperature anomalies in the Amazon, while its influence on NBP sensitivity to soil-moisture anomalies will be less pronounced.

4 Conclusions

In conclusion, in this study we assessed the main environmental drivers affecting the long-term and IAV projections of the Amazon basin carbon sink under the SSP5-8.5 high-radiative-forcing scenario with respect to 11 CMIP6 ESMs. Several important outcomes have been identified and confirmed.

Specifically, we recognized CO₂ fertilization as the predominant mechanism determining the long-term Amazon carbon sink trend and uncertainty under the range of atmospheric CO₂ concentrations of the SSP5-8.5 scenario (400–1135 ppmv). These results reaffirm, for the CMIP6 generation of ESMs, what was reported in a previous study by Huntzinger et al. (2017), i.e. the predominant role, under different research assumptions, of vegetation sensitivity to CO₂ in shaping the net carbon sink variability across CMIP5-generation land surface models (LSMs). Additionally, we disentangled the fundamental physical processes behind net carbon sink discrepancies across ESMs, highlighting dominant mechanisms affecting simulated carbon flux uncertainty for the Amazon basin ecosystem. Particularly, we show that the ensemble divergence of NBP under a future warming scenario is largely determined by GPP modulation caused by shortwave incoming radiation and uncertainty in the representation of heterotrophic respiration sensitivity to both soil moisture and temperature. Our multi-model ensemble approach expands the results obtained in previous research (Ma et al., 2021) by allowing for an explicit consideration of model uncertainty. By showing the strong uncertainties across the driving factors of heterotrophic respiration, we suggest that ESMs not only differ in the positive modulation of temperature acting on R_h (controlled by the Q₁₀ equation), but also largely disagree in the association between soil moisture and soil decomposition rates leading to respiration fluxes, as recently pointed out by Guenet et al. (2024).

Additionally, our results point towards a stronger ENSO-driven temperature impact on carbon sink anomalies for the vast majority of ESMs, compared to the effects associated with deficits in water availability. Accordingly, the CMIP6 multi-model ensemble considered shows a robust and statistically significant increase in the carbon sink sensitivity to ENSO-driven temperature anomalies under global warming in the region of the Amazon basin; as a consequence, climate change is likely to significantly diminish the Amazon ecosystem capacity to function as a carbon sink, further aggravating the atmospheric CO₂ burden. This outcome highlights the critical role that global-warming-induced changes in the dynamics of modes of interannual variability could have in the Amazon region, thus contributing to the expanding scientific literature on the topic (see for example Le, 2023; Yan et al., 2023). Overall, despite being able to, at least partially, understand and quantify the contributions of different driving factors in affecting Amazon carbon-cycle projections, important uncertainties still remain in the projections of the Amazon response to global warming in the CMIP6 generation of ESMs. On another note, we acknowledge some limitations affecting our study, particularly the application of a linear framework to describe interannual variability in carbon fluxes and ENSO-mediated impacts. A potential approach to overcome this limitation could be the adoption of machine-learning algorithms able to detect non-linearities and describe the complex spatio-temporal features of the Earth system (Reichstein et al., 2019). This, together with a more up-to-date and uniform representation of vegetation physiological and biogeochemical processes in the next generation of ESMs, will hopefully help generate more coherent and reliable projections of terrestrial carbon fluxes, helping us to understand the prospect of the Amazon basin under global warming scenarios.