Slow-down of the greening trend in natural vegetation with further rise in atmospheric CO2

Satellite data reveal widespread changes of Earth’s vegetation cover. Regions intensively attended to by humans are mostly greening due to land management. Natural vegetation, on the other hand, is exhibiting patterns of both greening and browning in all continents. Factors linked to anthropogenic carbon emissions, such as CO2 fertilization, climate change and consequent disturbances, such as fires and droughts, are hypothesized to be key drivers of changes in natural vegetation. A rigorous regional attribution at biome-level that can be scaled into a global picture of what is behind the observed changes 5 is currently lacking. Here we analyze different datasets of decades-long satellite observations of global leaf area index (LAI, 1981-2017) as well as other proxies of vegetation changes, and identify several clusters of significant long-term changes. Using process-based model simulations (Earth system and land surface models), we disentangle the effects of anthropogenic carbon emissions on LAI in a probabilistic setting applying Causal Counterfactual Theory. The analysis prominently indicates the effects of climate change on many biomes – warming in northern ecosystems (greening) and rainfall anomalies in tropical biomes 10 (browning). Our results provide only little support to previously published accounts of dominant global-scale effects of CO2


Introduction
Satellite observations reveal widespread changes in terrestrial vegetation across the entire globe. The greening and browning trends reflect changes in the abundance of green leaves, and thus, the rate and amount of photosynthesis. Plants modulate pivotal land-atmosphere interactions through the process of photosynthesis. Hence, changes in photosynthetic activity have immediate effects on the land-atmosphere exchange of energy (Forzieri et al., 2017), water (McPherson, 2007;Ukkola et al., 20 2016) and carbon (Poulter et al., 2014;Thomas et al., 2016;Winkler et al., 2019a, b). Several studies have reported that many biomes are largely greening, from Arctic tundra to subtropical drylands (Myneni et al., 1997;Nemani et al., 2003;Mao et al., 2016;Zhu et al., 2016;Chen et al., 2019;Winkler et al., 2019a). Others have identified regions of declining trends in leaf area (Goetz et al., 2005;Verbyla, 2011). The drivers underlying these long-term vegetation changes, however, remain under debate.
In the light of nearly forty years of continuous satellite observations, we reassess the driver attribution of natural vegetation 25 changes in a new cause-and-effect framework.
Anthropogenic vegetation, i.e. actively cultivated vegetation, and natural vegetation should be considered separately due to their distinct origins and properties. A recent study by Chen et al. (2019) reported that anthropogenic vegetation (35% of the global vegetated area) is greening due to human land management. The authors identified irrigation, multiple cropping, and the application of fertilizers and pesticides as the main drivers of leaf area enhancement (direct drivers). These results challenge 30 the conclusions of a previous study by Zhu et al. (2016) that attributed the global greening trend mostly to indirect drivers induced by CO 2 emissions, in particular, the CO 2 fertilization effect (70%).
Indirect drivers of vegetation changes usually include CO 2 fertilization and climatic change in the literature, both of which are consequences of rising atmospheric CO 2 concentration. The term "CO 2 fertilization" includes two effects of increased ambient CO 2 on the physiology of plants. First, elevated CO 2 in the interior of leaves stimulates carbon assimilation, which 35 enhances plant productivity and biomass (Leakey et al., 2009;Fatichi et al., 2016). Second, leaves adapt to the CO 2 -enriched atmosphere by lowering their stomatal conductance and, over time, stomatal density. As a consequence, water loss through transpiration decreases, resulting in increased water-use efficiency (ratio of carbon assimilation to transpiration rate; Ukkola et al., 2016;Fatichi et al., 2016). In theory, both effects should result in an expansion of leaf area, especially in environments where plant growth is constrained by water availability (Donohue et al., 2009(Donohue et al., , 2013Ukkola et al., 2016). 40 The radiative effect of CO 2 in the atmosphere induces climatic changes that can have both damaging or beneficial effects on the functioning of ecosystems. Temperature-limited biomes are expected to green (i.e. increase leaf area) due to warming and associated prolongation of the growing season (Park et al., 2016;Winkler et al., 2019a). But long-term drying (Zhou et al., 2014), as well as increased intensity and frequency of disturbances (Seidl et al., 2017) such as droughts (Bonal et al., 2016) and wildfires (Goetz et al., 2005;Verbyla, 2011), can induce regional vegetation browning trends. Regional greening and browning signals, is prone to overfitting, and assumes that linear correlation reflects causation (Hannart and Naveau, 2018). To overcome these limitations, we propose to use the Causal Counterfactual Theory which has recently been introduced to climate science (Pearl, 2009;Hannart et al., 2016;Hannart and Naveau, 2018). The method allows us to test if long-term greening/browning trends can be attributed to the effects of rising CO 2 in a probabilistic setting combining necessary and sufficient causation (Section 2.7). This is the first study that addresses vegetation browning as well as greening patterns across all major biomes, integrated 85 into a global picture. Greening is dominant in terms of areal fraction, but browning clusters are intensifying, primarily in the tropical forests that are biodiversity-rich and highly productive. We find that CO 2 fertilization is an important driver of greening in some biomes, but cannot be established as a dominant causal driver in many others. The strengthening browning trend identified in our study is most likely linked to climate changes, i.e. long-term drying and recurring droughts. Overall, our findings suggest that the emerging browning clusters in the highly productive ecosystems might be a precursor of a weakening 90 land carbon sink, which is not yet captured by the current land components of Earth system models.
The LAI3g datasets prior to 2000 were not evaluated due to a lack of required field data Chen et al., 2019).
After 2000, the quality of the LAI3g dataset was assessed through direct comparisons with ground measurements of LAI and indirectly with other satellite-data based LAI products, and also through statistical analysis with climatic variables such as 105 temperature and precipitation variability . Various studies used the predecessor LAI3gV0 and the related dataset of fraction of absorbed photosynthetically active radiation (FAPAR; Anav et al., 2013;Forkel et al., 2016;Zhu et al., 2016;Mao et al., 2016;Mahowald et al., 2016;Piao et al., 2014;Poulter et al., 2014;Keenan et al., 2016) and its successor LAI3gV1 (Winkler et al., 2019a, b;Chen et al., 2019).
Leaf area index is defined as the one-sided green leaf area per unit ground area in broadleaf canopies and as one-half the 110 green needle surface area in needleleaf canopies in both satellite observations and models (ESMs and LSMs). It is expressed in units of m 2 green leaf area per m 2 ground area. Missing values in the LAI3gV1 dataset are filled using the climatology of each 16-day composite during 1982-2017. We use the annual averaged LAI of each pixel in this study.
In addition to the GIMMS LAI3g product, we analyze the MODIS LAI record as well as four other long-term global remote sensing datasets: the Global Land Surface Satellites LAI product (GLASS LAI), the Global Mapping LAI product (GLOBMAP 115 LAI), the NDVI product from the Land Long Term Data Record (LTDR), and a new FAPAR product from the National Oceanic and Atmospheric Administration's (NOAA) Climate Data Record Program.
The MODIS LAI data analysed in this study are based on the combined Terra and Aqua MODIS LAI products (MOD15A2H and MYD15A2H) from Collection 6 (C6; Myneni et al., 2015a, b). These LAI datasets are provided at a 8-day temporal resolution with a 500 m sinusoidal projection covering the entire globe. After additional data quality assessment (for more 120 details, please see Chen et al., 2019), the two LAI datasets are aggregated into a 16-day composites by taking the mean of all valid LAI values. The data are then spatially aggregated to 1/20 degree spatial resolution and cover the period from 2000 to 2019.
The GLASS LAI dataset  is based on AVHRR, MODIS, and CYCLOPES reflectances and LAI products.
The full time series was generated using an artificial neural network (general regression neural networks) that has been trained 125 on the overlap period of AVHRR, MODIS and CYCLOPES reflectances and LAI products . We use the latest version of the GLASS LAI dataset, which covers the period from 1981 to 2018 and is provided at a 8-day temporal and a 1/20 degree spatial resolution.
The GLOBMAP LAI dataset (Liu et al., 2012) is a reconstruction of the historical AVHRR data by a quantitative fusion with MODIS data. The algorithm inverses a geometrical optical model to establish pixel-level relationships between AVHRR 130 and MODIS LAI for the overlapping period, which are then used to reconstruct AVHRR LAI back to the initial year of the record. We use the latest version of the GLOBMAP LAI dataset, which covers the period from 1981 to 2017 and is provided at a 15-day temporal and a 1/13.75 degree spatial resolution.
The NDVI dataset of NASA's LTDR (Pedelty et al., 2007) project is based on a reprocessing of long term AVHRR reflectances applying improved preprocessing techniques and atmospheric corrections used in the generation of MODIS datasets.

135
The preprocessing improvements include radiometric in-flight vicarious calibration for the visible and near-infrared channels and inverse navigation to relate an Earth location to each sensor instantaneous field of view (Pedelty et al., 2007). Atmospheric corrections include corrections for Rayleigh scattering, ozone, water vapor, and aerosols. We use the recently published version 5 (v5) of the LTDR NDVI dataset, which covers the period from 1981 to 2019 and is provided at a daily temporal and a 1/20 degree spatial resolution.

140
The FAPAR product from the National Oceanic and Atmospheric Administration's (NOAA) Climate Data Record Program provided by National Centers for Environmental Information (NCEI) is based on carefully calibrated and corrected land surface reflectances from AVHRR sensors (Claverie et al., 2016). The algorithm relies on artificial neural networks calibrated per different land cover types using the MODIS FAPAR dataset. We use the latest version of the FAPAR dataset from 1981 until 2019, which is provided at a daily temporal and a 1/20 degree spatial resolution. 145 We aggregate all data sets to annually averaged values and to spatially area-weighted averages for different biomes as defined by the mask (Section 2.2). The MODIS land cover classification does not contain the biome tundra, which is why we use in addition the land cover 160 product GLDAS2 / Noah version 3.3 that uses a modified IGBP classification scheme providing the classes Wooded, Mixed or Bare Ground Tundra (https://ldas.gsfc.nasa.gov/gldas/GLDASvegetation.php, hereafter GLDAS land cover) (Rodell et al., 2004). Accordingly, pixels originally of the classes Shrublands, Grasslands, Permanent Wetlands, or Barren, are converted to Tundra, if classified as Wooded, Mixed or Bare Ground Tundra in the GLDAS land cover product. The classes Woody Savannas and Savannas span vast areas across the globe in the MODIS land cover product. We use the GLDAS classification 165 for these pixels, but only for regions where the MODIS and GLDAS land cover products disagree. In doing so, we obtain a more accurate global land cover classification. Table S1 describes in detail how the fusion of the MODIS and GLDAS land cover products is realized.

Characterization of biomes & clusters of significant change
As a last step, we integrate the MODIS tree cover product MOD44B (MODIS/Terra Vegetation Continuous Fields Yearly L3 Global 250 m SIN Grid V006, https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod44b_v006) to 170 account for the underestimation of forested area in the MODIS land cover product. Areas with tree cover exceeding 10% are formally defined as forests (MacDicken et al., 2015). Thus, we set non-forest pixels in the MODIS land cover product above 10% tree cover to Boreal Forest in the high latitudes 50 • N/S. For tropical forest (25 • S -25 • N), we increase the threshold to 20% tree cover to allow for a realistic areal extent of savannas. The pixels in the bands 25 • N/S -50 • N/S remain unchanged, because the MODIS land cover product already realistically represents the forested area in these latitudes.

175
Table S1 provides a detailed overview on the conflation of MODIS land cover product, GLDAS land cover product and the MODIS Tree cover product. The final biome map (originally resolved at 0.05 • ) is regridded to the different resolutions of the AVHRR sensor and the models simulations (MPI-ESM and TRENDYv7) applying a largest area fraction remapping scheme.
Based on the observational LAI dataset we define various clusters for greening or browning in most biomes: North American Tundra (NAm Tundra), Eurasian Tundra (EA Tundra), North American Boreal Forests (NAm Brl F), Eurasian Boreal Forests Specifically, we used the initial CMIP6 release of the MPI-ESM version 1.2.01 (mpiesm-1.2.01-release, revision number 9234). The final CMIP6 version will include further bug fixes, which are expected to only slightly influence long-term sensitivities of simulated land surface processes. 205 We conducted historical simulations (all forcings) and three factorial experiments (all forcings except one): (a) all historical forcings except the physiological effect of CO 2 (No PE; increasing CO 2 does not affect the biogeochemical processes), (b) all historical forcings except the radiative effect of CO 2 (No RE; increasing CO 2 does not affect climate), and (c) all historical forcings except anthropogenic forcings (No CO 2 ). All experiments were preformed in ensemble-mode (6 realizations per experiment) using the latest CMIP6 forcing data . Individual realizations were initialized from different points in 210 time of a prolongation run of the official MPI-ESM1.2-LR pre-industrial control simulation. In doing so, we account for the influence of climatic modes (e.g. El Niño Southern Oscillation) as a source of uncertainty in simulating long-term changes.
The simulated time series were shifted by four years to maximize the overlap with the observational record of 1982-2017.

Land surface models: TRENDYv7
Land-surface models (LSMs) or dynamic global vegetation models (DGVMs) simulate key physical and biological key pro-215 cesses of the land system in interaction with the atmosphere. LSMs provide a deeper insight into the mechanisms controlling terrestrial energy, hydrological and carbon cycles, as well as the drivers of phenomena ranging from short-term anomalies to long-term changes (Sitch et al., 2015;Bastos et al., 2018). Here, we analyze the most recent TRENDY ensemble (version 7) comprising 13 state-of-the-art LSMs which vary in their representation of ecosystem processes. All models simulate vegetation growth and mortality, deforestation and regrowth, vegetation and soil carbon responses to increasing atmospheric CO 2 levels, 220 climate change and natural variability (Quéré et al., 2018). Some models simulate an explicit nitrogen cycle (allowing for potential nitrogen limitation) and account for atmospheric N deposition (Table A1 in Quéré et al., 2018). Most LSMs include the most important components of land-use and land-use changes, but they are far from representing all processes resulting from direct human land management (Table A1 in Quéré et al., 2018). A more detailed description of the TRENDYv7 ensemble, model-specific simulation setups and references can be found in Quéré et al. (2018, Table A4). 225 We use output from five simulations: all forcings (S3), physiological effect of CO 2 only (S1), radiative plus physiological effect of CO 2 (S2), land-use changes only (S4), and the control run (S0; no forcings: fixed CO 2 concentration of 276.59 ppm and fixed land-use map, loop of mean climate and variability from 1901-1920). The forcing data consist of observed atmospheric CO 2 concentrations, observed temporal patterns of temperature, precipitation, and incoming surface radiation from the CRU-JRA-55 reanalysis (Quéré et al., 2018;Harris et al., 2014), and human-induced land-cover changes and management 230 from an extensions of the most recent Land-Use Harmonization (LUH2) dataset (Hurtt et al., 2011;Quéré et al., 2018).
In this study, we only analyze output for the period 1982-2017 (matching the observational record) from models providing spatially gridded data for all five simulations. A few models provide LAI at the level of plant functional types (PFTs). We calculate the average value of all LAI values on PFT level multiplied by their land cover fraction for each grid cell. All model outputs were spatially regridded to a common resolution of 1 • based on a first-order conservative remapping scheme (Jones,235 1999).
The design of factorial simulations in TRENDYv7 and by the MPI-ESM are conceptually different. The MPI-ESM simulations were conducted using the counterfactual approach, i.e. all forcings are present except the driver of interest. TRENDYv7 provides simulations with different combinations of drivers as described above. To obtain comparability, we have to make the assumption that the absence of a specific driver has the same effect, in absolute values, as its sole presence. Thus, we process 240 the output of the simulations S1, S2, S3 and S4 to obtain the counterfactual setup as described above for MPI-ESM. This approach neglects possible synergy effects from simultaneously acting forcings. Also, it has to be noted that these simulations are only to some extent comparable between the two ensembles. For instance, in the MPI-ESM we can specifically determine the impact of the radiative effect of CO 2 , whereas TRENDYv7 uses observed atmospheric fields including changes induced from other drivers, such as non-CO 2 greenhouse gases.

245
For certain clusters, some models show unreasonable LAI changes and/or extreme inter-annual variability. To reduce the influence of these extreme models on the overall analysis, we apply a two-step filtering method for each cluster beforehand.
Models are excluded from the analysis, if they exceed three times the inter-annual variability of observations and/or show a drastic change (of either sign) of more than 250% between the start and end of the observational period. Further, we apply a weighting scheme based on the performance of the all-forcings run for each cluster. We calculate quartic weights based on 250 the distance between the simulated and observational estimate. These weights are applied when calculating the multi-model average and standard deviations for the factual and counterfactual runs.

Atmospheric CO 2 concentration
Global monthly means of atmospheric CO 2 concentration are taken from the GLOBALVIEW-CO2 product (for details see http://dx.doi.org/10.3334/OBSPACK/1002) provided by the National Oceanic and Atmospheric Administration/Earth System 255 Research Laboratory (NOAA/ESRL).

Processing of the gridded data
Areas of significant change in LAI are estimated using the non-parametric Mann-Kendall test, which detects monotonic trends in time series. In this study, we set the significance level to p ≤ 0.1. An alternative statistical test for trend detection (Cox-Stuart test;Sachs, 1997) yields approximately the same results. The trends are either calculated for time series on the pixel level or 260 for area-weighted large-scale aggregated time series (e.g. biome level).
We define greening (browning) either as a positive (negative) temporal trend, or for better comparison among models and observations as well as for a better global comparison across diverse biomes, we express these trends relative to the initial LAI level at the beginning of the observational record (average state from 1982-1984), denoted as Λ (% decade −1 ).
The calculation of yearly net changes in leaf area balances the effects from both statistically significant browning and 265 greening grid cells. For each cell, we multiply the estimated trends by the respective grid area. The net change is the sum of all grid cells, where areas of insignificant change are set to zero.
Models fairly accurately reproduce global patterns of vegetation greening, however, the fraction of browning is considerably underrepresented. Yet, we can only consider pixels with significant negative trends in LAI, in observations and models alike, and test models with respect to driver attribution of browning trends. Thus, the attribution of browning trends in this paper 270 exclusively refers to browning pixels only.
Models reveal biases in comparison to observations. To obtain informative results in the attribution analysis, we process the simulations to match the mean and variance of the observational time-series. Assuming additive and multiplicative biases in simulations, we apply the following corrections: This processing step does not affect the nature of simulated trends.

Causal Counterfactual Theory
The causal counterfactual approach is anchored in a formal theory of event causation developed in computer science (Pearl, 2009;Marotzke, 2019). Recently, a framework for driver attribution of long-term trends in the context of climate change has been introduced (Hannart et al., 2016;Hannart and Naveau, 2018), and increasingly gains popularity (Marotzke, 2019).

285
Through the use of this method we can ascertain the likelihood that a certain external forcing has caused an observed change in the Earth system. More precisely, we address the question of interest in a probabilistic setting, i.e. what is the probability that a given forcing (e.g. radiative effect of CO 2 ) has caused an observed long-term change in the system (e.g. greening of the Arctic).
In the following, we highlight the key ideas and relevant concepts of causal theory. A detailed description and formal 290 derivations can be found in (Pearl, 2009;Hannart et al., 2016;Hannart and Naveau, 2018). We define the cause event (C) as "presence of a given forcing" (i.e. the factual world that occurred) and the complementary event (C) as "absence of a given forcing" (i.e. the counterfactual world that would have existed in the absence of a given forcing; Hannart and Naveau, 2018). Further, we define the effect event (E) as the occurrence of a long-term change (here, greening or browning) and the complementary event (E) as the non-occurrence of a long-term change (i.e. no persistent vegetation changes). In making 295 use of numerical models, we can conduct factual runs comprising all forcings (i.e. historical simulations) as well as simulate counterfactual worlds by switching off a forcing of interest (i.e. all forcings except one). Based on an ensemble of simulations, either in a multi-model and/or multi-realizations setup, we derive the so-called factual (p 1 ) and counterfactual probability (p 0 ), which read p 1 = P {E|do(C)} and p 0 = P {E|do(C)}, respectively (Hannart and Naveau, 2018). More precisely, p 1 describes the probability of the event E in the real world where forcing C was present, whereas p 0 refers to the probability of the event 300 E in a hypothetical world where forcing C was absent. The notation do(·) means that an experimental intervention is applied to the system to obtain the probabilities (Hannart and Naveau, 2018).
The three distinct facets of causality can be established based on the probabilities p 1 and p 0 : PN refers to the probability of necessary causation, where the occurrence of E requires that of C but may also require other forcings. PS refers to the probability of sufficient causation, where the occurrence of C drives that of E but may not be required for E to occur. PNS describes the probability of necessary and sufficient causation, where PN and PS both hold (Hannart and Naveau, 2018). In other words, PNS may be considered as the probability that combines necessity and sufficiency. Thus, the 310 main goal is to establish a high PNS that reflects and communicates evidence for the existence of a causal relationship in a simple manner (Hannart and Naveau, 2018).
To obtain PNS, we follow the methodology described in detail in Hannart and Naveau (Hannart and Naveau, 2018) and derive cumulative distribution functions (CDF) for the factual and counterfactual worlds, denoted D 0 and D 1 , respectively.

315
Assuming a Gaussian distribution, PNS follows as where µ 1 and µ 0 refer to the mean response of all factual and all counterfactual runs, respectively. Σ denotes the overall uncertainty and is estimated based on all simulations, comprising factual, counterfactual, and centuries-long unforced (preindustrial) model runs (for details see Hannart and Naveau, 2018). Finally, the maximum of PNS determines the sought proba-320 bility of causation (Hannart and Naveau, 2018). We express probabilities using the terminology and framework defined by the IPCC (Mastrandrea et al., 2011;Hannart and Naveau, 2018). (i.e. including agriculturally dominated regions). The higher browning proportion in the extended record analyzed in this study indicates an intensification of leaf area loss in recent years. In the following, we take a closer look at different major biomes and their changes in LAI.

Earth's forests respond diversely throughout the satellite era
A global map of statistically significant trends in LAI (denoted Λ, Section 2.6) for natural vegetation reveals greening (Λ > 0) 335 and browning (Λ < 0) clusters across the globe (Figure 1). Temperate forests (Λ > 0: 56%) and Eurasian boreal forests (Λ > 0: 53%) exhibit extensive regions of increasing LAI, and thereby, contribute the largest fraction to the enhancement of leaf area on the planet ( Table 2)  interestingly, in that they show changes of reversed sign; Figure 1).

355
Leaf area loss occurs primarily in densely vegetated biomes (i.e. forests), which outweighs leaf area gain in rather sparsely vegetated regions (i.e. grasslands). For instance, vigorously greening areas of circumpolar tundra result in a leaf area gain of 8.74 × 10 3 km 2 yr −1 , which is almost outbalanced fourfold by a leaf area loss of 34.31 × 10 3 km 2 yr −1 in the browning regions of the tropical forests (Table 2). To assess the responses of different biomes to rising CO 2 in more detail, we iteratively calculate statistically significant LAI trends for different time windows with advancing initial year (i.e. 1982, 1983, ..., 2000),  Table S3). Critically, the tropical forests display the sharpest transition from a substantial net gain of 24.11 × 10 3 km 2 yr −1 (Table 2) to a comparably strong net loss of leaf area (-18.42 × 10 3 km 2 yr −1 ; Table S3). To address the temporal development of positive and negative changes in leaf area in more detail, we calculate time series of area-weighted averages of LAI ( Figure 2a). We find that browning of natural vegetation occurs at a considerably higher level of LAI (on average ∼1.85) than greening (on average ∼1.32). Throughout the observational period, these two time series of opposite trends converge 385 towards a LAI of 1.6 ( Figure 2a). This convergence of greening and browning is not only evident in terms of their LAI level ( Figure 2a), but also in their proportions (inset in Figure 1). The time series of anthropogenic vegetation on the other hand, aggregated for positive and negative Λ separately, are both confined to a comparable low LAI level (on average between 1 and 1.25). We next investigate the global LAI distributions of negative and positive changes and their development over time.
Comparing distributions of the earlier (1982)(1983)(1984) with those of the more recent years (2015-2017) reveals that browning 390 primarily occurs at a high (5-6) and a medium level of LAI (1-2.5; Figure 2b). Greening, however, is occurring almost entirely at low levels of LAI between 0-1.5. As a consequence, the global area-weighted averages of the browning and greening regions are approaching one another (dashed versus solid vertical lines in Figure 2b), as also depicted by the time series (Figure 2a).
Overall, these results suggest a homogenization of Earth's natural vegetation in terms of LAI texture with rising CO 2 . This homogenization becomes prominent when we compare the distributions of negative and positive Λ over time using a Q-Q 395 plot (quantile-quantile; Figure 2c). The relationship between the quantiles is skewed to the left at higher LAI (positive Λ on x-axis, negative Λ on y-axis), because browning is prevalent in high LAI regions. Over time, the quantiles of the greening and browning distributions are approaching the 1-1 line (representing identical distributions), emphasizing their convergence.

The majority of models reproduce the observed convergence of greening and browning trends
Thus far, we have described the diverse long-term changes of natural vegetation across all continents and throughout the satellite 400 era. We next investigate the underlying mechanisms driving these greening and browning trends and use the fully-coupled MPI-ESM and the TRENDYv7 ensemble of observation-driven LSMs (Section 2.3 and 2.4). First, we ask if these models capture the observed behavior of natural vegetation under rising CO 2 . MPI-ESM reproduces the observed browning of high LAI and the greening of low LAI regions, however, the levels of LAI do not match the observations ( Figure S2). Historical simulations of TRENDYv7 (here 13 models) also show pronounced changes in vegetation, but exhibit a diverse behavior among the models 405 (results not shown for brevity). Seven LSMs reproduce observed converging trends of greening and browning, whereas the other six models show divergent trends. All TRENDYv7 models are driven with identical atmospheric forcing fields, hence, these six models most likely lack or incorrectly represent key processes of ecosystem functioning. In general, simulated greening patterns are comparable to observations (Murray-Tortarolo et al., 2013;Sitch et al., 2015;Mahowald et al., 2016), but browning, especially in the North American boreal forests, is underestimated (Sitch et al., 2015). We use the framework of Counterfactual Causal Theory to attribute changes in LAI to a given driver in a probabilistic setting (Pearl, 2009;Hannart et al., 2016;Hannart and Naveau, 2018). Note that the causal relationships in this approach are determined based on the predictions of the models for the all-forcings (also referred to as factual) and factorial runs (also referred to as counterfactual), and therefore the causality results may reflect biases or misrepresentations in the models. Based 420 on the factual and counterfactual runs, we derive Probabilities of causation that combines Necessity and Sufficiency of each factor (PNS; see Section 2.7 for details). When aggregated to area-weighted global averages (i.e. Earth greening trend), the observed estimate (∼ 1.08 % decade −1 ) and the factual MPI-ESM estimate (∼ 1.14 % decade −1 ) are comparable, whereas the multi-model average of the TRENDYv7 ensemble is an overestimate (∼ 1.79 % decade −1 ; Figure 3c). Omitting CO 2 -induced climate change (no radiative effect of CO 2 , No RE) does not have a strong effect in the MPI-ESM (∼ 1.04 % decade −1 ), i.e.

425
the estimate does not differ considerably from the factual run. The TRENDYv7 models indicate that the positive trend in LAI can be explained by climate change to some extent (∼ 1.21 % decade −1 ). However, the PNS values for the radiative effect of CO 2 are generally rather low (Fig. 3d), implying that the probability of the radiative effect of CO 2 acting as a sufficient and necessary causal driver of the globally aggregated LAI trend signal is rather low. The opposite is the case, when the physiological effect of CO 2 (No PE) is excluded. Both model setups agree that almost no positive trend in LAI is present in a 430 world without the CO 2 fertilization effect (MPI-ESM: ∼ 0.18 % decade −1 , TRENDYv7: ∼ 0.08 % decade −1 ; both estimates are lower than internal variability of ∼ 0.49 % decade −1 ). Note that the term variability here refers to a broader concept of variability that includes several components, such as climate variability, inter-model variability, and, if applicable, variability in observations, adapted from the approach introduced by Hannart and Naveau (2018, see also Section 2.7 for details).
As a consequence, high PNS can be established: The physiological effect of CO 2 has in the case of MPI-ESM likely (68%) 435 and in the case of TRENDYv7 very likely (91% ) caused the positive trend of global LAI in recent decades (Figure 3d). This result is in line with Zhu et al. (2016) who reported that 70% of global greening is attributable to CO 2 fertilization. Removing both effects of CO 2 results in slight negative trends, probably due to land use practices (e.g. deforestation; Figure 3c).

The global signal switches to a minor negative trend in the second half of the observational period
Natural vegetation shows a slight negative trend for the period 2000-2017 (∼ -0.4 % decade −1 ; Figure 3e). This estimate is 440 within the range of internal variability, and thus, should be interpreted with caution. Note, that the net change in leaf area is still positive when considering only significantly changing pixels (inset in Figure 1). To provide confidence in this result, we analyze three additional remote sensing datasets for LAI (MODIS-LAI, GLASS-LAI and GLOBMAP-LAI) as well as for the normalized difference vegetation index (LTDR-NDVI) and for the fraction of absorbed photosynthetic active radiation (NCEI-FAPAR), both proxies for leaf area changes (see Section 2.1 for details). We calculate time series of changes relative to the  Wang et al. (2020) suggests that global CO 2 fertilization has declined in recent years by analyzing various observational datasets and highlights that land-surface models are not reproducing this development, mainly due to the under-representation of nutrient limitation. While these results are consistent with ours, we are not convinced that one can infer a decline or saturation of the CO 2 fertilization effect from these observational datasets. Rather, we argue that countervailing effects associated with the radiative effect of increasing CO 2

465
(climatic changes, e.g., increase in atmospheric dryness and changes in water availability), as discussed below in more detail, become more pronounced and increasingly reduce vegetation productivity.
Overall, driver attribution at the global scale, as described above, and also in Zhu et al. (2016), neglects the heterogeneity of natural vegetation and the possibility that divergent responses of different natural biomes might cancel out. To account for this omission, we identify eleven clusters of significant change and derive probabilities of causation for each driver across different vegetation types ( Figure 5).

Temperate forests prosper with rising CO 2 while tropical forests are increasingly under stress
Forests in temperate climates exhibit a strong positive trend in LAI (∼ 2.53 % decade −1 ), which is also seen in the models, albeit slightly overestimated (MPI-ESM: ∼ 3.18 % decade −1 , TRENDYv7: ∼ 2.69 % decade −1 ; Figure S3). The physiological effect of CO 2 is the main driver with high PNS (85% for MPI-ESM, 80% for TRENDYv7; Figure 5). The trends are slightly 475 weaker when only analyzing the second half of the observational period, but the overall result does not change. Observed warming might have additionally contributed to enhanced vegetation growth (e.g. growing season extension; Piao et al., 2011;Park et al., 2016), however, it is not identified as an important driver by models. Most temperate forests are in developed countries, and thus, have been managed in a sustainable manner for several decades (Currie and Bergen, 2008). It is conceivable, that some of the positive trends in LAI could be attributed to forest management or regrowing forests (Pugh et al., 2019), 480 however, this is not captured by the models (i.e. trends are negative when complete CO 2 forcing is absent; Figure S3).
The response of tropical forests to rising CO 2 is more complex. The signal over the entire observational period is slightly positive (∼ 0.3 % decade −1 ), however, it is within the range of internal variability. Therefore, no robust driver attribution is possible ( Figure 5 and Figure S4). TRENDYv7 models show strongly opposing responses of LAI to the different effects of CO 2 : LAI decreases when the physiological effect is omitted, but increases when the radiative effect is omitted. MPI-

485
ESM shows qualitatively the same responses, but less pronounced ( Figure S4). For the second half of the satellite record, the observed trend switches sign to a strong negative trend (∼ -1.4 % decade −1 ). The models reproduce this tendency, but the multi-model average of the TRENDYv7 ensemble is still positive. During the same time period, the opposing reactions to CO 2 in the factorial runs are more strongly marked ( Figure S4). These results suggest that browning caused by CO 2 -induced climate change is compensated by greening affiliated with the CO 2 fertilization effect at the biome level. Based on these findings, we 490 hypothesize that the physiological effect of CO 2 is strong in models and outbalances the negative effect of climate change in the tropical forests (Kolby Smith et al., 2016). As a consequence, the all-forcings simulations fail to reproduce the observed patterns of strengthening vegetation browning in the tropics (Zhou et al., 2014;Song et al., 2018). Because of the demonstrated limited predictive power of the models in simulating the vegetation response to climatic changes, we also rely more heavily on the published literature in the following discussion of our results. in sea surface temperatures in the tropical North Atlantic (Marengo et al., 2008(Marengo et al., , 2011Xu et al., 2011). Whereas the ENSO-driven droughts peak in northern hemispheric winter, thus during the wet season, the non-ENSO droughts happened during the dry season (July -September), when tropical ecosystems are more vulnerable to negative rainfall anomalies.
These intense and frequent droughts have diverse impacts on tropical ecosystems (Bonal et al., 2016), the most prominent being an increase in wildfires and tree mortality. Recently, perennial legacy effects have been identified which lead to persistent 505 biomass loss in the aftermath of severe droughts Yang et al., 2018). For instance, some regions were still recovering from the impact of the megadrought of 2005 when the next major drought began in 2010 . Maeda et al. (2015) found that these extreme events are also capable of disrupting hydrological mechanisms, which can lead to long-lasting changes in the structure of Amazonian ecosystems. Such droughts and associated wildfires are predicted to increase in frequency (Cai et al., 2014) and intensity (Fasullo et al., 2018) as a consequence of the ENSO-related amplification 510 of heat waves, but also due to the projected warming of the tropical North Atlantic (Munday and Washington, 2019).
In addition to these episodic disturbances, long-term changes in climate also affected the tropical forests in the Amazon region. Rising surface air temperatures have considerably increased atmospheric water vapor pressure deficit (VPD), which has a negative effect on vegetation growth (Yuan et al., 2019). Moreover, we find that precipitation has steadily decreased during the dry season (July -September, Figure S5 and S6) based on the latest version of the ECMWF reanalysis for the last 515 forty years (ERA5; Dee et al., 2011). This rainfall deficit and the identified lengthening of the dry season (Fu et al., 2013) exacerbate vegetation water stress during dry seasons and favor conditions for wildfires. The slight increasing trend in wet season precipitation (February -April) most likely cannot compensate for the water loss and its impact during the dry season ( Figure S5). Overall, the intensification of the dry season and the recurring droughts cause long-term browning trends , in line with our results of intensified browning of Amazonian forests ( Figure S6). 520

Drying trend in central African humid forests
African tropical forests have been experiencing a long-term drying trend since the 1970s (Malhi and Wright, 2004;Asefi-Najafabady and Saatchi, 2013;Zhou et al., 2014). In contrast to South America, the steady decline in rainfall is seen during both dry and wet seasons ( Figure S5). The origin of this decreasing trend in year-round rainfall is still under debate. Precipitation in equatorial Africa is expected to increase under climate change (Weber et al., 2018), so it is hypothesized that this trend is 525 associated with the Atlantic Multidecadal Oscillation and/or changes in the West African Monsoon system (Asefi-Najafabady and . Long-term drying in rainforests could also be connected to the physiological effect of rising CO 2 . Recently, it has been demonstrated that the reduction in stomatal conductance and transpiration induces a drier, warmer, and deeper boundary layer, resulting in a decline in local rainfall (Langenbrunner et al., 2019). Regardless of what the causes may be, this long-term water deficiency most likely has led to the most pronounced cluster of vegetation browning in Earth's tropical 530 forests (∼ 174 × 10 3 km 2 net loss of leaf area in the time period of 2000-2017). No robust attribution is possible with the set of models analyzed in this study, since they fail to capture this substantial decrease in leaf area in the all forcing runs ( Figure S7).
In the case of the TRENDYv7 models, this finding is particularly noteworthy as they are driven with observed precipitation changes: The spatial patterns of negative trends in LAI and dry season precipitation in the Central African tropical forests coincide to a large extent ( Figure S5).
Interestingly, the MODIS record does not exhibit this browning cluster , though it has been reported in other independent observational datasets (Zhou et al., 2014). Also, atmospheric CO 2 inversion studies have identified negative trends in carbon uptake for this region (Fernández-Martínez et al., 2019), which corroborates our results based on the LAI3g dataset.

Tropical forests in Oceania are afflicted by deforestation
Although we exclude direct anthropogenic land-cover changes ( Figure S1, Table S1) as well as abrupt changes (Mann-Kendall test for monotonic trends, Section 2.6), the LAI trend maps nevertheless show characteristic deforestation patterns, e.g. the so-called "arc of deforestation" in the Amazon region ( Figure S6; Aldrich et al., 2012). Hence, deforestation practices may explain some part of the observed gradual browning of the Amazon (Song et al., 2015) and African tropical forests (Mayaux et al., 2013;Tyukavina et al., 2018).

545
In Oceania, however, deforestation appears to be a crucial driver of the observed browning in the pristine tropical forests.
Significant negative trends align strongly with patterns of drastic deforestation during recent decades, described in detail by Stibig et al. (2014, in comparison to Figure 1). As opposed to Central Africa and the Amazon region, climate changes are unlikely to be the key driver of browning regions in Oceania. There, precipitation, although highly variable in the dry season, appears to increase ( Figure S5) and the increase in VPD is rather minor in tropical forests (Yuan et al., 2019). and TRENDYv7: ∼ 2.08 % decade −1 ), which can mostly be attributed to amplified warming of the temperature-limited northern high latitudes (PNS = 71% for TRENDYv7, PNS = 44% for MPI-ESM; Figure S8). North American boreal forests exhibit a negative response to the effects of rising CO 2 , which has amplified over the last two decades (∼ -0.95 % decade −1 , 2000-555 2017). Models do not reproduce the dominant browning pattern ( Figure S9), which is most likely connected to inadequate representation of disturbances (Sitch et al., 2015). Several studies have proposed that browning has occurred as consequence of droughts, wildfire, and insect outbreaks in the North American boreal forests (Goetz et al., 2005;Sitch et al., 2015;Beck and Goetz, 2011;Kurz et al., 2008). Macias Fauria and Johnson (2008) showed that the frequency of wildfires is strongly related to the dynamics of large-scale climatic patterns (Pacific Decadal Oscillation, El Niño Southern Oscillation, and Arctic Oscillation) 560 and thus, cannot be tied conclusively to anthropogenic climate change. However, there is also evidence that the residing tree species suffer from drought stress induced by higher evaporative demand as the temperature rises (Verbyla, 2011). Moreover, models lack a representation of the asymmetry in tree species distribution between North America and Eurasia, which could explain their divergent reactions to changes in key environmental variables (Abis and Brovkin, 2017). Further observational evidence for the browning of North American boreal forests and the associated decline in net ecosystem productivity can also 565 be inferred from CO 2 inversion products (Fernández-Martínez et al., 2019;Bastos et al., 2019).
Tundra ecosystems also reveal a dipole-type development between North America and Eurasia, however with a reversed sign. Hence, North American tundra is strongly greening (Observations: ∼ 4.23 % decade −1 , MPI-ESM: ∼ 4 % decade −1 , and TRENDYv7: ∼ 4.51 % decade −1 ), which is virtually certain (PNS = 99% for TRENDYv7) and about likely as not (PNS = 51% for MPI-ESM) caused by warming ( Figure S10). The trend decreases for the period 2000-2017, which could be linked 570 to the warming hiatus in the years 1998-2012 (Bhatt et al., 2013;Ballantyne et al., 2017;Hedemann et al., 2017). This is in line with the observed slow down in tundra greening due to short-term cooling after volcanic eruptions (Lucht et al., 2002).
Eurasian tundra show a positive trend for the years 1982-2017, but a reversal in trend sign for the years 2000-2017 ( Figure   S11). Models exhibit some evidence of a strengthening browning signal, but fail to capture the full extent of the emerging browning clusters seen in observations. If we only consider the grid cells that show significant browning in observations and 575 models, we are able to conduct a robust driver attribution. According to the TRENDYv7 ensemble, the browning cluster in Eurasian tundra can very likely be attributed to CO 2 induced climate change (PNS = 93%, PNS = 47% for MPI-ESM). These results are in line with studies showing that tundra ecosystems are susceptible to warm spells during growing season (Phoenix and Bjerke, 2016) and to frequent droughts (Beck and Goetz, 2011). The asymmetry between Eurasia and North America can be explained by changes in large-scale atmospheric circulation. Eurasia is cooling through increased summer cloud cover,

580
whereas North America is warming through more cloudless skies (Bhatt et al., 2013(Bhatt et al., , 2014. Also linkages between regional Arctic sea ice retreat, subsequent increasing ice-free waters, and regional Arctic vegetation dynamics have been postulated (Bhatt et al., 2014).

Vegetation in arid climates is greening, except in South America
Non-forested greening clusters beyond the high northern latitudes coincide with semi-arid to arid climates (Park et al., 2018).

585
The Northern Sub-Saharan African savannas and grasslands greened extensively in recent decades (∼ 4.63 % decade −1 ; Figure   S12), which is reproduced by the observation-driven TRENDYv7 models (∼ 4.55 % decade −1 ), and is likely caused by climatic changes (PNS = 68%). No robust attribution is feasible based on the MPI-ESM simulations. However, it is noteworthy, that the fully-coupled Earth system model points to climate change as having a negative effect in these regions, thus, not reproducing the observed increase in rainfall ( Figure S12). This provides evidence for the hypothesis that African precipitation anomalies 590 are not induced by rising CO 2 , but rather follow a multidecadal internal climatic mode (Asefi-Najafabady and . Internal variability in LAI changes is strong in the Southern African grasslands and savannas, and thus, no robust long-term change can be identified ( Figure S13). It has been shown that shrublands in the more southern regions are greening in response to increased rainfall (Fensholt and Rasmussen, 2011). In general, the literature suggests that greening and browning patterns in arid climates are mainly driven by precipitation anomalies (Fensholt and Rasmussen, 2011;Fensholt et al., 2012;Gu et al., 595 2016; Adler et al., 2017). Close resemblance arises when comparing the spatial patterns of precipitation trends throughout the satellite era (Adler et al., 2017) with significant changes in vegetation in arid environments, especially so in the African continent. Decreased rainfall in arid South America coincides with strong browning clusters (Fensholt et al., 2012). This is in disagreement with the expected strong manifestation of CO 2 fertilization in water-limited environments (Ukkola et al., 2016).
Australian Shrublands show a persistent positive LAI trend (∼ 3.84 % decade −1 ), intermittently perturbed by climatic 600 extreme events (e.g. strong anomalous rainfall with subsequent extensive vegetation greening in 2011, Figure S14; Poulter et al., 2014). Models reproduce the steady greening of Australia, but no robust driver attribution is feasible due to strong internal variability. However, both model setups point to the physiological effect of CO 2 as the dominant driver ( Figure S14).
These results are congruent with recent studies (Donohue et al., 2009;Ukkola et al., 2016) that show CO 2 fertilization enhanced vegetation growth by lowering the water limitation threshold.

605
Grasslands in the cool arid climates exhibit persistent positive trends (∼ 2.03 % decade −1 , Figure S15). Simulated estimates are in the range of the observations (MPI-ESM: ∼ 2.33 % decade −1 and TRENDYv7: ∼ 1.81 % decade −1 ). Our analysis suggests that the positive response of cool arid grasslands to rising CO 2 can be explained by the physiological effect of CO 2 (PNS = 85% for TRENDYv7, PNS = 88% for MPI-ESM). These ecosystems are dominated by C 3 -type plants (Still et al., 2003), which are susceptible to CO 2 fertilization (Sage et al., 2012), thus, consistent with our results. In the warm arid areas, 610 C 4 -type grasses dominate (Still et al., 2003), which are less sensitive to the physiological effects of CO 2 (Sage et al., 2012).
As discussed above, vegetation changes there are mostly driven by precipitation anomalies, although CO 2 fertilization might also contribute to a limited extent (Sage et al., 2012).  T e m p e r a t e F o r e s t s E u r a s i a n B o r e a l F o r e s t s N -A f r i c a n S a v a n n a s -G r a s s l a n d s A u s t r a l i a n S h r u b l a n d s C o o l G r a s s l a n d s N -A m e r i c a n T u n d r a N -A m e r i c a n B o r e a l F o r e s t s E u r a s i a n T u n d r a S -A f r i c a n S a v a n n a s -G r a s s l a n d s In this paper we examine nearly four decades of global LAI changes under rising atmospheric CO 2 concentration. We find that 615 the Earth's greening trend is weakening and clusters of browning are beginning to emerge, and importantly, expanding during the last two decades. With one exception, all analyzed satellite observation datasets confirm these results, but with different signal strengths. Leaf area is primarily decreasing in the pan-tropical green belt of dense vegetation. Leaf area gain is occurring mostly in sparsely vegetated regions in cold and/or arid climatic zones, and in temperate forests. Thus, vegetation greening is occurring mainly in regions of low LAI, whereas browning is seen primarily in regions of high LAI. Consequently, these 620 opposing trends are decreasing the texture of leaf area distribution in natural vegetation.
We identify clusters of greening and browning spread across all continents and conduct a regional, i.e. biome-specific, driver attribution based on factorial model simulations. The results suggest that the physiological effect of CO 2 (i.e. CO 2 fertilization) is the dominant driver of increasing leaf area only in temperate forests, cool arid grasslands and likely the Australian shrublands.
A cause-and-effect relationship between CO 2 fertilization and greening of other biomes could not be established. This finding 625 questions the study by Zhu et al. (2016) that identified CO 2 fertilization as the globally prevailing driver of the Earth's greening trend. We find that many clusters of greening and browning bear the signature of climatic changes. The greening of Sub-Saharan grasslands and savannas is consistent with increase in rainfall. Climatic changes, primarily warming and drying, determine the patterns of vegetation changes in the northern ecosystems, i.e. greening of Eurasian boreal forests and North American tundra, but also emerging browning trend in the Eurasian tundra. Models fail to capture the browning of North American boreal forests.

630
Models suggest rising CO 2 has compensatory effects on LAI in the tropical forests. Climatic changes induce browning, which is opposed by greening due to a strong physiological effect in the models. Hence, if the physiological effect of CO 2 is "turnedoff", models simulate the emerging browning trend in the tropics comparable to observations. Our analysis of changes in rainfall during the satellite age underpins climate changes as the main cause of tropical forest browning: recurrent droughts and decline in dry season precipitation in the Amazon as well as long-term drying trends in Africa.

635
Models represent a simplified view of the real world reduced to its essential processes. Some of these processes are underrepresented or lacking in the current generation of land-surface models. Whether they are driven with observed climatic conditions or operate in a fully-coupled Earth system model, they fail to capture the full extent of adverse effects of rising CO 2 in natural ecosystems. In particular, the deficiency of reproducing the observed leaf area loss in North American boreal and in pantropical forests -biomes which account for a large part of the photosynthetic carbon fixation -has considerable implications 640 for future climate projections. Thus, it is important to focus model development not only on a better representation of disturbances such as droughts and wildfires, but also on revising the implementation of processes associated with the physiological effect of CO 2 , which currently offsets browning induced by climatic changes.
Another vital issue for future research is the impact of large-scale climatic anomalies on vegetation. All three major clusters of browning are hypothesized to be associated with temperature or precipitation anomalies modulated by climatic modes. Many 645 droughts in the Amazon were attributed to El Niño events (Bonal et al., 2016). The long-term drying trend in tropical Africa is possibly connected to the Atlantic Multidecadal Oscillation (Asefi-Najafabady and . Likewise, disturbances in North American boreal forests are likely controlled by an interplay between large-scale climatic patterns (Pacific Decadal Oscillation, El Niño Southern Oscillation, and Arctic Oscillation;Macias Fauria and Johnson, 2008). Little is known about how these large-scale patterns might change in a warming climate. Current Earth system models struggle to simulate these 650 climatic modes and related precipitation patterns, which is likely rooted in their coarse spatial resolution. New tools, such as high-resolution simulations or large ensembles, offer possibilities to study these phenomena.
Overall, our study suggests that the Earth largely greened in the 1980s and 1990s as rising CO 2 triggered mainly LAIincreasing effects, e.g., by warming the high northern latitudes and overall more carbon allocation through CO 2 fertilization.
However, as CO 2 continues to rise, the system appears to be entering or has entered a regime in which LAI-decreasing effects 655 are amplified, i.e., climatic changes associated with rising atmospheric CO 2 concentration become more pronounced and have stronger adverse effects in various ecosystems. In addition, plant sensitivity to CO 2 fertilization may already be saturating, as recently suggested by Wang et al. (2020), but this aspect remains controversial.
We show that the effects of rising CO 2 on LAI are not comparable across the biomes, as are the impacts on the ecosystems.
Regarding biodiversity, the consequences of leaf area loss in tropical forests that harbor the most diverse flora and fauna of the 660 planet are not compensated for by leaf area gain in temperate and arctic ecosystems. A similar caveat is in order with respect to the carbon cycle, e.g. an additional leaf in the tundra does not offset the reduction in primary productivity of a leaf lost in the tropical rainforest. Thus, our results indicating loss of tropical leaf area should be of concern. A recent study suggested that the tropical forests have already switched to being a net source of carbon, also considering land-use emissions (Baccini et al., 2017). The uncertainty in future projections is large, ranging from a stable CO 2 fertilization-driven carbon sink to a collapse 665 of the system at a certain CO 2 concentration (Cox et al., 2000). Concerning leaf area, the models project a steady greening of the tropical forests in the high-end CO 2 emissions scenario (business-as-usual) and a slight browning in the low-end scenario (mitigation) by the end of the century . Altogether, the tropical forests have the potential to crucially influence the evolution of climate throughout the 21 st century and should be a vital issue for future research. ences provided in respective Methods section. Processed data and analysis scripts are available from the corresponding author upon request and will also be published in public repositories together with this article.