Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification — TEI-tagged view

A. L. BurrellpersName, J. P. EvanspersName, M. G. De KauwepersName (2020date). “Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification.” Nature Communicationstitle·j.

Object of study — desertificationterm·dom of global drylandsterm·dom and its drivers, 1982-2015

Research question — how much anthropogenic climate changeterm·dom and land useterm·dom each drive desertificationterm·dom

Contribution — the first observation-based attributionterm·key of desertificationterm·dom including CO2 fertilization and land-use break points

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Anthropogenic climate changeterm·dom has driven over 5 million km2 of drylandsterm·dom towards desertificationterm·dom Abstract Drylandsterm·dom cover 41% of the earth’s land surface and include 45% of the world’s agricultural land. These regions are among the most vulnerable ecosystems to anthropogenic climate and land useterm·dom change and are under threat of desertificationterm·dom. Understanding the roles of anthropogenic climate changeterm·dom, which includes the CO 2 fertilizationterm·dom effect, and land useterm·dom in driving desertificationterm·dom is essential for effective policy responses but remains poorly quantified with methodological differences resulting in large variations in attributionterm·key. Here, we perform the first observation-based attributionterm·key study of desertificationterm·dom that accounts for climate change, climate variabilityterm·dom, CO 2 fertilizationterm·dom as well as both the gradual and rapid ecosystem changes caused by land useterm·dom. We found that, between 1982date and 2015date, 6% of the world’s drylandsterm·dom underwent desertificationterm·dom driven by unsustainable land useterm·dom practices compounded by anthropogenic climate changeterm·dom. Despite an average global greening, anthropogenic climate changeterm·dom has degraded 12.6% (5.43 million km 2 ) of drylandsterm·dom, contributing to desertificationterm·dom and affecting 213 million people, 93% of who live in developing economies.

Introduction Land degradationterm·key is a systemic global problem 1 – 4 but the scale of the problem is disputed, with global estimates of degraded areas ranging from <10 to >60 million km 2 5 . Changes in vegetation in drylandsterm·dom are predominantly caused by two factors: (i) anthropogenic climate changeterm·dom (ACCterm·dom), which includes both changes in water availability driven by trends in precipitation and increases in temperature 6 , 7 , as well as increased water use efficiency (carbon gain per unit of water lost) in response to rising atmospheric CO 2 8 ; and (ii) land useterm·dom (LUterm·dom) practices, including grazing, cropping and deforestation 2 , 9 . Unsustainable LUterm·dom is considered the primary negative driver of drylandterm·dom degradationterm·dom 9 – 11 . The impact of climate change (CC) on drylandsterm·dom is also generally thought to be negative, with some studies suggesting that anthropogenic forcing has already increased arid areas 12 – 14 .

Despite evidence for LU-induced degradationterm·dom and the studies that find increased aridification over drylandsterm·dom, satellite estimates of vegetation greennessterm·dom (a proxy for net primary productivityterm·dom (NPPterm·dom)) show a significant global increase since 1980 10 . The key drivers of this global increase in apparent vegetation productivity are the vegetation’s response to rising CO 2 8 , 15 , increases in rainfall and temperature 16 , 17 and LUterm·dom 10 . Model simulations which prescribe LUterm·dom, attribute almost all of the trend in satellite-derived greening to CO 2 fertilizationterm·dom 15 , while satellite-derived models that do not account for CO 2 , explicitly find either climate or LUterm·dom as the dominate factor 10 , 17 . Neither approach explicitly accounts for rapid ecosystem change (break points) in their proportioning of the relative contributions of each driver. This can lead them to miss or underestimate rapid changes driven by processes like extreme fires, deforestation, reforestation, changes in agricultural policy, etc. 18 – 21 . Disentangling the roles of climate (temperature and precipitation), CO 2 and LUterm·dom thus remains a key challenge 22 and has been identified as a key knowledge gap by the United Nations Convention to Combat Desertificationterm·dom 2 (UNCCD), the Intergovernmental Panel on Climate Change 23 (IPCC), and the Intergovernmental Science–Policy Platform on Biodiversity and Ecosystem Services 3 .

Here we quantified the scale of global desertificationterm·dom, which both the UNCCD and IPCC define as degradationterm·dom in arid, semi-arid, and dry sub-humid areas 23 . These definitions further define degradationterm·dom as the long-term reduction or loss of biological productivity among other things. Here we identify areas undergoing long-term reductions in vegetation in drylandterm·dom areas, hence the desertificationterm·dom according to relevant international conventions, using the satellite-based GIMMSv3.1g Normalized Difference Vegetation Index (NDVI) data. We calculated the overall vegetation change using a non-parametric trend analysis applied to peak growing season NDVI (NDVI max ). We then attributed this change to CO 2 , climate variabilityterm·dom (CVterm·dom), CC, and LUterm·dom using a modified version of the Time Series Segmented Residual Trendsterm·mth (TSS-RESTRENDterm·mth) method 19 , 24 . This approach quantifies the effect of interannual CVterm·dom as well as long-term changes in climate and CO 2 fertilizationterm·dom in addition to ecosystem break points caused by LUterm·dom (see “Methods”). To quantify uncertainties, we used a 12-member ensemble made up of statistical model runs performed using a combination of observation-based gridded datasets (four precipitation and three temperature datasets). We show that 6% of drylandterm·dom areas have undergone desertificationterm·dom since 1982date with a further 20% of drylandterm·dom areas being at high risk of future desertificationterm·dom as a result of unsustainable LUterm·dom practices and ACCterm·dom.

Results and discussion The extent and drivers of drylandterm·dom vegetation change Globally, of the 44.5 million km 2 of drylandsterm·dom, 6% of these areas experienced desertificationterm·dom (i.e., significant negative change in NDVI max ), 41% showed significant greening (i.e., significant positive change), and 53% had no significant change between 1982date and 2015date (Fig. 1a ). The mean (±1 SD) of the area-weighted drylandterm·dom vegetation change, as represented by the change in NDVI max was 0.031 ± 0.053. We estimated the scale of desertificationterm·dom to be 2.70 million km 2 , which is significantly below a previous estimate of ~10.5 million km 2 over the same region, but over a different time window (1982date and 2003) 1 . A large part of this discrepancy can be attributed to climatic differences in the end dates of the studies (2003 vs. 2015date), with increased rainfall over regions including the Sahel and India 25 , 26 . This large difference between our estimate and this existing drylandterm·dom degradationterm·dom estimate highlights that time-series vegetation trend analysis is sensitive to the start and end conditions 17 , 18 . For this reason, understanding what is driving the observed vegetation change is more important than the current directions of vegetation change for projecting future changes and vulnerabilities. It should also be noted that although the amount of land we estimate to have experienced desertificationterm·dom has decreased by ~70% between these two estimates, the number of people impacted has only decreased by ~25% (250 million to 189 million) 27 . Fig. 1 Primary drivers of vegetation changes between 1982date and 2015date. a The observed drylandterm·dom vegetation change from 1982date to 2015date measured using the change in p Normalized Difference Vegetation Index (ΔNDVI max ). Non-dryland and hyper-arid regions are masked in dark gray and areas where the change is insignificant ( α FDR = 0.10) or smaller than the error in the sensors (±0.001) are masked in white. b Map showing the largest absolute attributed driver (CO 2 , land useterm·dom, climate change, and climate variabilityterm·dom). Non-dryland regions are masked in dark gray and the magnitude of other factors is not considered. c The mean and magnitude (mean absolute value) of the different drivers of change in NDVI max between 1982date and 2015date. The error bars show the SD of the area-weight grid cells. d Frequency distribution function for the change in NDVI max (1982date–2015date) attributed to land useterm·dom, which shows that there are opposing changes at local scales which cancel out when averaged globally.

Figure 1b shows that globally, CO 2 fertilizationterm·dom was the largest absolute attributed driver of drylandterm·dom vegetation change in 44.1% of areas, followed by LUterm·dom (28.2%), CVterm·dom (14.6%), and then CC (13.1%). However, when averaged globally (Fig. 1c ), the per-pixel contribution of CO 2 (0.021 ± 0.011) was much larger than the contribution from CVterm·dom (0.006 ± 0.020), CC (−0.002 ± 0.023), or LUterm·dom (0.005 ± 0.032). The relative contribution (67.8% CO 2 , −5.6% climate, 15.5% LUterm·dom) fall within the range of global estimates calculated using a model based factor analysis (70.1 ± 29.4% CO 2 fertilizationterm·dom, 8.1 ± 20.6% climate, 3.7 ± 14.7% LUterm·dom) 15 , despite a difference in the study domains. It should be noted here that model based factor analysis did not quantify the role of CVterm·dom 15 , which we find accounts for 19.4% of the observed global drylandterm·dom greening between 1982date and 2015date.

If only the global mean effect size is considered, climate and LUterm·dom seem to have a very small impact compared with CO 2 fertilizationterm·dom, which seemingly contradicts well-documented evidence of LUterm·dom and climate impacts 9 , 11 , 16 , 17 , 28 , and the spatial patterns shown in Fig. 1b . For example, a recent satellite-based study, which did not consider the role of CO 2 , attributed 60% of observed global land changes to LUterm·dom activities and the remaining 40% to other factors including climate 10 . For comparison, we repeated our ensemble analysis, removing the role of CO 2 fertilizationterm·dom (Supplementary Fig. 1 ), and found that 60.4% of global drylandterm·dom vegetation change would have been attributed to LUterm·dom and 39.6% to CC and variability combined, which is consistent despite a difference in the study domains and the attributionterm·key methods (see Supplementary Text 1 ). This result underlines the need to explicitly consider the positive role of CO 2 as a driving mechanism for change in drylandterm·dom ecosystems 8 . When CO 2 is included in the analysis, the mean effect of LUterm·dom and climate are small, because there are roughly equal areas of positive and negative change that largely cancel out when averaged globally (Fig. 1c, d ), a result which holds even if we assume a different level of vegetation response to elevated CO 2 (Supplementary Fig. 2 ). This means it is also important to consider the magnitude of the different drivers (Fig. 1c ).

The impact of anthropogenic climate changeterm·dom Combining the change due to CO 2 fertilizationterm·dom and CC provides a quantification of the role of ACCterm·dom in recent desertificationterm·dom (Fig. 2a ). Globally, we found that ACCterm·dom had a positive (greening effect) over the study period (NDVI max : 0.019 ± 0.027). Although broadly positive, ACCterm·dom also had a desertifying effect across 12.55% (5.43 million km 2 ) of drylandsterm·dom areas. Hotspots where ACCterm·dom had a desertifying effect include parts of the western United States, eastern Brazil, Iraq, Syria, Jordan, Kazakhstan, Uzbekistan, Mongolia, and Australia. Crucially, the negative effects of ACCterm·dom are disproportionately felt by poorer nations with 85% of the 213.4 million people impacted living in developing or newly industrialized countries. However, a negative ACCterm·dom forcing does not guarantee an area experienced desertificationterm·dom (Fig. 3a ). Only 13.8% (0.75 million km 2 ) of areas with a negative ACCterm·dom forcing, experienced significant desertificationterm·dom ( α FDR = 0.10) and in only 2.27% (0.015 million km 2 ) of the areas experiencing desertificationterm·dom, did we find that climate was the sole negative driver. Fig. 2 The contribution of anthropogenic climate changeterm·dom to vegetation change from 1982date to 2015date. a The contribution of anthropogenic climate changeterm·dom (Climate Change + O 2 ) component to the change in vegetation between 1982date and 2015date (ΔNDVI max ). Non-dryland and hyper-arid regions are masked in dark gray and areas where the change is insignificant ( α FDR = 0.10) or smaller than the error in the sensors (±0.001) are masked in white. b The mean area-weighted pixel effect size of the drivers of observed vegetation broken down by observed vegetation change direction. Positive indicates greening and negative is desertificationterm·dom. Error bars show the SD. c Frequency distribution function for the change in NDVI max (1982date–2015date) attributed to by anthropogenic climate changeterm·dom. Fig. 3 Desertificationterm·dom risk and regional drivers. a Regions experiencing, or at risk of experiencing, desertificationterm·dom. Areas with a significant negative change in vegetation ( α FDR = 0.10) are classified “Desertificationterm·dom,” areas where the anthropogenic climate changeterm·dom (ACCterm·dom: CO 2 + Climate Change) and land-use (LUterm·dom) components are both negative but the change in vegetation is not significant are classified as “LUterm·dom and ACCterm·dom.” Areas where the anthropogenic climate changeterm·dom has had a negative effect are classified as “ACCterm·dom” and areas where land useterm·dom had a negative impact classified as “LUterm·dom.” b map of the regional subdivisions used. c The mean per-pixel ΔNDVI max and its drivers broken down by regions shown in b . The error bars show the SD.

Drivers of desertificationterm·dom In the 2.70 million km 2 of drylandsterm·dom that experienced desertificationterm·dom, a negative LUterm·dom component was the primary driver in 79.9% and a contributing factor across 99.0% of areas (Fig. 2b ). Even though the average impact of CC (NDVI max : −0.004 ± 0.030) and CVterm·dom (−0.002 ± 0.024) are much smaller than LUterm·dom (NDVI max : −0.040 ± 0.034), climate remains an important driver of desertificationterm·dom. Ecosystems that are experiencing reduced water availability or drought conditions are much more vulnerable to degradationterm·dom from LUterm·dom and vice versa, with the negative effects compounding 7 . For example, over parts of Central Asia we observed negative changes in both CC and LUterm·dom (Fig. 4 ), which is consistent with the strong evidence of long-term degradationterm·dom driven by unsustainable LUterm·dom practices resulting in the well-documented Aral Sea disaster 26 . Similarly, the negative impacts of decreased rainfall over the semi-arid Caatinga forest of Brazil has amplified the effects of widespread deforestation and grazing intensification. The decrease in rainfall in South America results from both CC (mean precipitation anomaly 1982date–2015date ≈ −0.2) as well as a negative phase of CVterm·dom (mean precipitation anomaly ≤ 0 and mean temperature anomaly ≥ 0) from 2009 to 2014 (Supplementary Figs. 3 and 4 ). For further discussion of these regions and comparison with regional studies, see Supplementary Text 2.1 . In addition to the drylandsterm·dom experiencing desertificationterm·dom, there are additional 12.0 million km 2 and 507 million people, living in areas where the desertifying effect of LUterm·dom has been offset by a positive ACCterm·dom signal (Fig. 3 ). These regions, along with the 7.2% of areas with a negative CC, but no significant vegetation change, are at the highest risk of future desertificationterm·dom. Fig. 4 The drivers of global vegetation change. The changes in NDVI max between 1982date and 2015date (ΔNDVI max ) attributed to ( a ) CO 2 fertilizationterm·dom, ( b ) climate variabilityterm·dom, ( c ) climate change, and ( d ) land useterm·dom. Non-dryland and hyper-arid regions are masked in dark gray. Areas where the change did not meet the multi-run ensemble significance criteria detailed in the “ Methods,” or are smaller than the error in the sensors (±0.001) are masked in white.

Drivers of drylandterm·dom greening We also observed widespread global greening, with 18.0 million km 2 of drylandsterm·dom having a significant positive vegetation change (Fig. 1a ). CO 2 was the largest driver of this change (Figs. 1b and 2b ), in line with previous findings of a dominant CO 2 fertilizationterm·dom effect on global vegetation greening 8 , 15 . Where our results differ from previous findings is in highlighting the importance of LUterm·dom and CVterm·dom. Unlike model based approaches which prescribe LUterm·dom 15 and have large variability in the simulated response of vegetation to climate 29 , our approach empirically determines the impact of climate and LUterm·dom on a per-pixel basis using observations. In the regions that experienced greening, we find that CO 2 was the largest attributed driver in ~40% of areas followed closely by LUterm·dom (~38%), CVterm·dom (~13%), and CC (~8%). The importance of CVterm·dom and LUterm·dom is especially apparent when considering regional drivers, with one or both playing a large role in the observed greening in the Sahel, India, China and Australia (Figs. 3c and 4 , and Supplementary Discussion 2.2 ).

We used an ensemble approach to minimize the uncertainty caused by observational datasets and structural change detection to account for the ecosystem break points driven by processes like deforestation observed in ~20% of areas (Supplementary Fig. 5 ). Our estimate of the total attributable change (CO 2 + LUterm·dom + CC + CVterm·dom) varied from the observed vegetation change by only 3% and when mapped spatially, the observed vegetation change and the total attributable change show very consistent patterns of greening and browning (Supplementary Fig. 6 ). Furthermore, our greening attributionterm·key results are consistent with regional studies done in the Sahel, India, China, and Australia (Supplementary Discussion 2.2 ).

In summary, understanding the causes of drylandterm·dom degradationterm·dom is an important and necessary step in targeting mitigating action that can reduce the impact of CC and prevent widespread desertificationterm·dom. Our change detection and attributionterm·key approach highlights the importance of accounting for the role of CO 2 and accurately quantifying the impact of LUterm·dom when considering potential drivers of change in drylandterm·dom ecosystems. Our results show that, despite widespread vegetation greening, 6% of areas that have undergone desertificationterm·dom mostly over western Asia and South America. This desertificationterm·dom directly effects 190 million people. In addition, we showed that unsustainable LUterm·dom practices or ACCterm·dom has placed 20% of drylandsterm·dom at high risk of desertificationterm·dom. This impacts 580 million people with the risk experienced disproportionately by low socioeconomic countries. Overall, our results highlight the importance of understanding what is driving the vegetation change for projecting impacts, because without this understanding there is a high risk that mitigation strategies will fail to prevent desertificationterm·dom.

Methods Quantifying desertificationterm·dom There is no universally agreed upon definition of desertificationterm·dom 5 , 23 , 30 . Here we use the UNCCD definition of land degradationterm·key, which is a reduction or loss of the biological or economic productivity resulting from various factors, including climatic variations and human activities, with desertificationterm·dom being any land degradationterm·key in drylandterm·dom ecosystems 2 . Drylandsterm·dom are defined by the UNCCD to be areas with an Aridity index < 0.05 or >0.65 2 . Historically, this has been measured using a linear trend applied to a satellite-derived vegetation proxy 1 . In this study, we used the growing season maximum NNDVI (NDVI max ) as a proxy of vegetation growth. For most regions, peak growing season NDVI was determined using the maximum value in a calendar year. However, pixels where peak the occurred in December, the January and February NDVI values of the subsequent year are considered part of the previous year’s growing season.

NDVI max has been found to have a highly significant correlation with NPPterm·dom in a large range of different drylandterm·dom ecosystem 31 – 33 . The Desertificationterm·dom chapter of the 2019 IPCC report on CC and LUterm·dom both, the UNCCD definition of desertificationterm·dom, and NDVI max as a proxy of vegetation growth 23 . It should be noted that the UNCCD definition of desertificationterm·dom and the trend in vegetation data used to measure it will not identify processes such as shrub encroachment, over intensification of agriculture, or the invasions by non-native species, which have been linked to degradationterm·dom but can cause increases in proxies such as NDVI 23 , 30 , 34 .

To produce a comparable estimate of desertificationterm·dom, we used a per-pixel non-parametric trend method (Theil–Sen slope estimator and Spearman’s ρ significance test), applied to the satellite-derived GIMMSv3.1g NDVI max dataset 35 . The global drylandterm·dom vegetation change (Obs) was calculated as the difference between the expected values (E) at the start (1982date) and the end (2015date) of the time series ( Obs = E 2015date − E 1982date ) . It should be noted that this is identical to multiplying the annual trend by the length of the time series, in this case 34 years. We report all variables using the difference between expected values at the start and end of the time series (ΔNDVI max ) rather than an as an annual trend to account for the ecosystem break points that are detected in some locations in the middle of the time series.

This study used version 3.1 of the 1/12° Global Inventory for Mapping and Modeling Studies (GIMMS) 35 NDVI dataset, which spans 1982date–2015date. Although the shorter temporal but higher spatial resolution datasets from newer sensors such Moderate Resolution Imaging Spectroradiometer (MODIS) do offer advantages, the shorter temporal record poses a serious issue in drylandterm·dom ecosystems. The natural variability in drylandterm·dom ecosystems is greatly impacted by decadal climate modes the most significant of which is El Niño Southern Oscillation (ENSO) 36 , 37 . The intensity of ENSO events varies significantly and since 1980 there have been three extreme El Niño (1982date, 1997, and 2015date) events that significantly impacted drylandterm·dom regions around the world 38 . Even with its almost 20-year record, MODIS has only captured one of these events (2015date) compared with the three present in GIMMS record. It is for this reason that GIMMS remains the most widely used dataset for vegetation trend detection and attributionterm·key studies 15 .

Accounting for the CO 2 fertilizationterm·dom effect To attribute the change in NDVI to the CO 2 effect on plant productivity between 1982date and 2015date, we used a theoretical relationship that links the increase in photosynthesis to increasing CO 2 39 (Eq. 1 ). 1 GPP ( rel ) ≈ c a − Γ * c a 0 + 2 Γ * c a + 2 Γ * c a 0 − Γ * where GPP (rel) is the relative CO 2 assimilation rate (%), c a is the atmospheric CO 2 concentration (µmol mol −1 ), and Γ* is the CO 2 compensation point in the absence of dark respiration (µmol mol −1 ). We set c a 0 to the CO 2 concentration in 1980 40 (~339 µmol mol −1 ) and Γ* = 40 (µmol mol −1 ).

Franks et al. 39 argued that the longer term response of plants to increasing CO 2 follows the ribulose 1,5-bisphosphate regeneration-limited rate (see also McMurtrie et al. 41 ). Accordingly, this relationship implies a conservative response to CO 2 (as plants may actually follow the Rubisco-limited rate when calculated on a intercellular CO 2 concentration ( C i ) basis during the period of 1982date–2015date 42 ) and ignores any indirect effects (i.e., increased water availability due to stomatal closure, which may extend the growth period in drylandsterm·dom, or interactions with seasonal rainfall 43 ). This approach has previously been advocated as a plausible assumption to correctly estimate gross primary productivity (GPP) using satellite light-use efficiency models 44 . We then assume that there has been no change in ratio of GPP to autotrophic respiration (Ra) during this period (1982date–2015date) and, as a result, the relative change in GPP equates to the relative change in NPPterm·dom based on Eq. 1 . During the period 1982date–2015date, global air temperatures have risen, which may have led to an increase in Ra 45 . However, increasing temperature has also increased carbon uptake, and both GPP and Ra have been shown to acclimate to the prevailing temperatures 46 , 47 , meaning that it does not necessarily follow that the GPP : Ra ratio has changed.

We apply the nonlinear CO 2 relationship (Eq. 1 ) to the raw NDVI data (NDVI obs ) to produce a scaled NDVI estimate (NDVI adj ) that excludes the CO 2 fertilizationterm·dom effect using Eq. 2 to relate the relative change in NPPterm·dom to a relative change in NDVI. 2 NPPterm·dom obs NPPterm·dom base ≈ NDVI obs NDVI adj where NPPterm·dom obs is the NPPterm·dom at the observed atmospheric CO 2 concentration ( c a ), NPPterm·dom base is the NPPterm·dom given the same climate conditions but an atmospheric CO 2 concentration of c a 0 , NDVI obs is measured NDVI value, and NDVI adj NPPterm·dom given the same climate conditions but an atmospheric CO 2 concentration of c a 0 . Equation 2 was used to calculate a NDVI adj value for every in the full NDVI obs time series with the atmospheric CO 2 concentrations taken from the IPCC historical forcing data 40 . This approach assumes that NPPterm·dom and NDVI are linearly related: .. where b ≈ 0 and m varies spatially. The linear relationship between NPPterm·dom and NDVI has been observed with both field estimates of NPPterm·dom 31 , 32 , 48 and estimates of NPPterm·dom derived from remote-sensing platforms 1 , 49 , 50 . However, this assumption of linearity breaks down in densely vegetated regions where NDVI saturates and in biomes with very low above-ground biomass, where the spectral characteristics of the bare soil influences NDVI values 51 , 52 . As we have excluded hyper-arid and non-water-limited ecosystems, which is a standard practice for studies on desertificationterm·dom (for more information, see ref. 23 ), we expect this assumption of linearity to be robust for our analysis (areas with Aridity index < 0.05 or >0.65 are masked from analysis).

The change in vegetation (ΔNDVI max ) attributed to the CO 2 fertilizationterm·dom effect was calculated by first taking the difference between peak growing NDVI with and without CO 2 fertilizationterm·dom effect (NDVI obs − NDVI adj ). Similar to the calculation of the observed ΔNDVI max , the non-parametric Theil–Sen slope estimator and Spearman’s ρ test for significance was then applied to these values. A time-series plot of the mean global NDVI obs and NDVI adj is shown in Supplementary Fig. 7 to highlight temporal nature of this attributionterm·key.

Determining the impact of climate and land useterm·dom After the NDVI was scaled to remove the CO 2 effect using the relationship from Franks et al. 39 , a statistical approach was used to attribute the change in NDVI adj to climate (both CVterm·dom and CC) and LUterm·dom. We used the recently developed TSS-RESTRENDterm·mth method, which allows vegetation changes due to LUterm·dom to be separated from those driven by CC and variability 19 , 24 , 53 .

Drylandterm·dom ecosystems have large natural interannual CVterm·dom. To separate the impact of climate and LUterm·dom, previous studies have generally fitted a statistical relationship between climate and vegetation, then used the trends in that relationship or its residuals to quantify LUterm·dom impacts 54 , 55 . However, LUterm·dom can have both a gradual impact on vegetation through processes such as grazing, which are captured by these methodologies 54 , 56 , and abrupt impact through processes such as deforestation, which cause these methods to break down 18 .

TSS-RESTRENDterm·mth differs from existing drylandterm·dom trend attributionterm·key methods in that it is able to capture both the long-term trends and the step changes in NDVI that occur in regions where ecosystems have experienced significant structural changes 19 . To do this TSS-RESTRENDterm·mth incorporates a phenological change detection method 57 to identify structural changes in the ecosystem, which manifest as break points in the NDVI time series. We used TSS-RESTRENDterm·mth v2.15, which has been updated to use both precipitation and temperature 24 , to calculate the Vegetation Climate Relationship (VCR) using the per-pixel optimal precipitation and temperature accumulation periods 19 , 24 . TSS-RESTRENDterm·mth was applied to the NDVI adj , with the LUterm·dom driven component calculated using an ordinary least squared regression between the residuals of the VCR and time, accounting for any detected structural changes to the ecosystem. A similar approach was presented by the IPCC 23 , although that approach does not separate CC from CVterm·dom and also assumes that all drylandterm·dom plants follow a C3 photosynthetic pathway.

Separating climate change and climate variabilityterm·dom The TSS-RESTRENDterm·mth method separates the effects of LUterm·dom from the combined effects of climate (variability and change) using VCR. To separate the effect of CC and CVterm·dom, the observed climatology (the per-pixel accumulated precipitation and temperature data) was calculated for the period 1962 to 2015date. A 20-year leading edge smoothing window was then applied to this observed climatology to remove the interannual CVterm·dom. The long-term trend caused by CC was determined using the Theil–Sen slope estimator 58 applied to the smoothed data and the results were used to detrend the observed climatology.

Using the per-pixel VCR, the total climate driven NDVI (NDVI CL ) and the NDVI due to CVterm·dom (NDVI cvterm·dom ) were calculated using the observed climatology and detrended climatology, respectively. The difference between NDVI CL and NDVI CVterm·dom is the change in NDVI max attributable to CC (NDVI CC ). The non-parametric Theil–Sen slope estimator and Spearman’s ρ test for significance was then applied to NDVI CVterm·dom and NDVI CC to get the change attributable to CVterm·dom and CC, respectively. The influence of Other Factors (OF), which could not be modeled, was calculated using OF = Obs − ( CO 2 + LUterm·dom + CVterm·dom + CC ) .

When discussing regional drivers of vegetation change in the main text and in Supplementary Discussion 2 , we report the mean climate anomaly rather than the accumulated precipitation and temperature values to allow comparison of different accumulation and offset periods of different pixels. For the observed accumulated precipitation and temperature, the anomaly was calculated on a per-pixel basis using: 3 z n = x n − μ obs. σ obs. where z = anomaly, n is the year, x = observed value, μ = mean of the per-pixel accumulated precipitation or temperature, σ = SD of the per-pixel accumulated precipitation or temperature. The temperature and precipitation anomaly attributed to CC was calculated using: 4 z n = β × ( n − n 0 ) σ obs. where β is the per-pixel trend in accumulated precipitation or temperature, n 0 is the first year of the analysis (1982date). The temperature and precipitation anomaly attributed to CVterm·dom was calculated using: 5 z n = x (adj) n − μ obs. σ obs. where x (adj) is the detrended precipitation or temperature. Regional breakdowns of the time series of observed, CC- and CV-driven precipitation and temperature anomaly are shown in Supplementary Figs. 3 and 4 , respectively.

Calculating the impact of Anthropogenic Climate Changeterm·dom In this study, we use the term ACCterm·dom to refer to both changes in water availability driven by trends in precipitation and increases in temperature 14 , 23 , as well as increased water use efficiency (carbon gain per unit of water lost) in response to rising atmospheric CO 2 8 , 15 . The impact of ACCterm·dom) is calculated using ACCterm·dom = CO 2 + CC. Although we acknowledge that this considering CO 2 fertilizationterm·dom and CC together is not common in the literature, these two things are aspects of the same anthropogenic cause; hence, we have reason to discuss them together. It should be noted that although the methodology used in this study is able to capture the effects of other anthropogenic greenhouse gases like methane in the CC component, we cannot separately quantify any direct impact that changes in the amount of these gasses will have on vegetation.

Accounting for dataset uncertainties and different photosynthetic pathways (C3 vs. C4) BurrellpersName et al. 53 showed that using an ensemble of TSS-RESTRENDterm·mth runs using different climate datasets improves the accuracy and minimizes the impact of errors and biases in the individual datasets. Climate data are relatively poorly sampled in drylandterm·dom regions, which can amplify the documented discrepancies between different datasets 59 , 60 We used two 12-member ensembles with matched runs made using TSS-RESTRENDterm·mth analysis performed using every combination of four precipitation and three temperature datasets (see Table 1 for details). The first ensemble assumes all plants are C3 and respond to elevated CO 2 , and the second ensemble assumes all plants are C4 with no eCO 2 response (see Supplementary Figs. 1 and 2 ). Table 1 Table of gridded datasets. Dataset Description References Vegetation Global Inventory for Mapping and Modeling Studies (GIMMSv3.1 g) NDVI, 15 day at 1/16° aggregated to monthly using the max of the valid values 35 Precipitation University of East Anglia Climate Research Unit TS v. 4.01 (CRU4 p ) Precipitation, monthly at 0.5° 72 – 76 The Climate Hazards group Infrared Precipitation with Stations v2.0 (CHIRPS) Precipitation, monthly at 0.05° 77 , 78 Multi-Source Weighted-Ensemble Precipitation (MSWEP) Precipitation, daily at 0.25° 79 TerraClimate Precipitation, monthly at 1/24° 80 Temperature University of East Anglia Climate Research Unit TS v. 4.01 (CRU4 T ) Temperature, monthly at 0.5° 66 – 70 TerraClimate Temperature, monthly at 1/24° 80 National Oceanic and Atmospheric Administration (NOAA) Climate Prediction Center (CPC) Temperature, daily at 0.5° 81 Additional datasets TerraClimate Potential Evapotranspiration (PET) used for the calculation of Aridity index (P/PET) 80 Global change in net primary productivityterm·dom (1981-2003), data from the Food and Agriculture Organization (FAO) Change in NPPterm·dom from 1982-2003 1 Gridded Population of the World version 4 (GPWv4) Gridded Population data 27 United Nations Development Programme Human Development Index (HDI) National HDI 82 North American Carbon Program (NACP) Global C3 and C4 SYNergetic land cover MAP (SYNMAP) C3/C4 vegetation fraction 61 A third 12-member ensemble, which accounts for the relative fraction of C3 and C4 plants, were calculated by taking the weighted mean of matched runs in ensemble one and two where the weights were the per-pixel fractions of C3 and C4 plants. Estimates of the relative fraction of C4 present in each pixel were derived from the matching 0.5° pixel in the North American Carbon Program (NACP) Global C3 and C4 SYNergetic land cover MAP (SYNMAP) 61 . The p -values for each ensemble member were combined using the same weights and the Stouffer’s Z -score method. All the results presented in the main paper are from this C4 adjusted ensemble.

We make the assumption that in drylandterm·dom ecosystems, precipitation controls the amount of foliage cover. As a result, increases in water use efficiency with increasing CO 2 will lead to increases in foliage cover in plants that use the C3 photosynthetic pathway 8 while plants that use the C4 photosynthetic pathway are not expected to show this response 62 . Our assumption that C3 plant respond strongly to increased atmospheric CO 2 , and that C4 plants do not respond, which is consistent with theory and short-term CO 2 enrichment experiments 62 . However, recent surprising results from a long-term CO 2 manipulation experiment in a drylandterm·dom ecosystem have shown much higher responses in C4 plants 62 . For this reason, the results of both the C3 and C4 12-member ensembles are included in the Supplementary Material for those interested parties (see Supplementary Figs. 1 , 2 , 8 , and 9 ).

Determining statistical significance For both Obs trend and CO 2 -driven change in NDVI, the Spearman’s rank correlation co-efficient test was used to measure statistical significance for each pixel 63 . In order to determine field significance and account for the multiple testing problem the Benjamini–Hochberg procedure was then applied to these p-values to control the false discovery rate (FDR) (α FDR = 0.10) 64 . As for the uncertainty in the approach we used to measure the CO 2 fertilizationterm·dom effect, the C3 and C4 12-member ensembles included in the Supplementary Material provide an estimate of the upper and lower bounds of CO 2 responses (see Supplementary Figs. 1 and 2 ).

For LUterm·dom, CVterm·dom, and CC there are 12 members in each ensemble. The p -values of these members were combined using the Fisher’s combined probability test and then the Benjamini–Hochberg procedure was applied to these p -values to determine statistical significance ( α FDR = 0.10). In addition, we also applied the IPCC protocol for determining ensemble significance and agreement 65 . For a pixel to be significant under this protocol, >50% of ensemble members must find a significant change ( α FDR = 0.10) and, of those significant runs, >80% must agree on the direction of change. If a pixel fails either the overall significance test or the IPCC protocol, the estimate of that component is masked for that pixel. Similar criteria have also been applied to ensemble breakpoint detection (>50% of runs must find a significant breakpoint, 80% of which must be in a three-year window) the results of which are included in supplementary material.

All the Climate datasets where interpolated from their native resolution to the 1/12° grid of the GIMMSv3.1 g datasets using the First-Order Conservative Remapping 66 in CDO 67 . Additional datasets where used to aid the interpretation and discussion our results. All dataset where converted to the GIMMS grid to allow for per-pixel comparison. To separate CC from CVterm·dom, a 20-year leading edge moving window was used where the value for a given year is the mean of the previous 20 years. For this reason, it was necessary to have climate data that goes back to 1960, which not all datasets do. Those climate datasets with insufficient temporal record were extended using TERRACLIMATE precipitation and temperature, and the Delta Bias Correction method 68 .

Supplementary information Supplementary Information Peer Review File