# Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification

Jul 31, 2020

### Quantifying desertification

There is no universally agreed upon definition of desertification5,23,30. Here we use the UNCCD definition of land degradation, which is a reduction or loss of the biological or economic productivity resulting from various factors, including climatic variations and human activities, with desertification being any land degradation in dryland ecosystems2. Drylands are defined by the UNCCD to be areas with an Aridity index < 0.05 or >0.652. Historically, this has been measured using a linear trend applied to a satellite-derived vegetation proxy1. In this study, we used the growing season maximum NNDVI (NDVImax) 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.

NDVImax has been found to have a highly significant correlation with NPP in a large range of different dryland ecosystem31,32,33. The Desertification chapter of the 2019 IPCC report on CC and LU both, the UNCCD definition of desertification, and NDVImax as a proxy of vegetation growth23. It should be noted that the UNCCD definition of desertification 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 degradation but can cause increases in proxies such as NDVI23,30,34.

To produce a comparable estimate of desertification, 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 NDVImax dataset35. The global dryland vegetation change (Obs) was calculated as the difference between the expected values (E) at the start (1982) and the end (2015) of the time series (({mathrm{{Obs}}} = E_{2015} – E_{1982})). 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 (ΔNDVImax) 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 1982–2015. 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 dryland ecosystems. The natural variability in dryland 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 (1982, 1997, and 2015) events that significantly impacted dryland regions around the world38. Even with its almost 20-year record, MODIS has only captured one of these events (2015) 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 attribution studies15.

### Accounting for the CO2 fertilization effect

To attribute the change in NDVI to the CO2 effect on plant productivity between 1982 and 2015, we used a theoretical relationship that links the increase in photosynthesis to increasing CO239 (Eq. 1).

$${mathrm{{GPP}}}_{({mathrm{{rel}}})} approx left[ {frac{{left( {c_a – {Gamma}^ ast } right)left( {c_{a0} + 2{Gamma}^ ast } right)}}{{left( {c_a + 2{Gamma}^ ast } right)left( {c_{a0} – {Gamma}^ ast } right)}}} right]$$

(1)

where GPP(rel) is the relative CO2 assimilation rate (%), ca is the atmospheric CO2 concentration (µmol mol−1), and Γ* is the CO2 compensation point in the absence of dark respiration (µmol mol−1). We set ca0 to the CO2 concentration in 198040 (~339 µmol mol−1) and Γ* = 40 (µmol mol−1).

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

We apply the nonlinear CO2 relationship (Eq. 1) to the raw NDVI data (NDVIobs) to produce a scaled NDVI estimate (NDVIadj) that excludes the CO2 fertilization effect using Eq. 2 to relate the relative change in NPP to a relative change in NDVI.

$$frac{{mathrm{{NPP}}}_{{mathrm{{obs}}}}}{{mathrm{{NPP}}}_{mathrm{{base}}}} approx frac{{mathrm{{NDVI}}}_{mathrm{{obs}}}}{{mathrm{{NDVI}}}_{mathrm{{adj}}}}$$

(2)

where NPPobs is the NPP at the observed atmospheric CO2 concentration (ca), NPPbase is the NPP given the same climate conditions but an atmospheric CO2 concentration of ca0, NDVIobs is measured NDVI value, and NDVIadj NPP given the same climate conditions but an atmospheric CO2 concentration of ca0. Equation 2 was used to calculate a NDVIadj value for every in the full NDVIobs time series with the atmospheric CO2 concentrations taken from the IPCC historical forcing data40. This approach assumes that NPP and NDVI are linearly related: .. where b ≈ 0 and m varies spatially. The linear relationship between NPP and NDVI has been observed with both field estimates of NPP31,32,48 and estimates of NPP derived from remote-sensing platforms1,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 values51,52. As we have excluded hyper-arid and non-water-limited ecosystems, which is a standard practice for studies on desertification (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 (ΔNDVImax) attributed to the CO2 fertilization effect was calculated by first taking the difference between peak growing NDVI with and without CO2 fertilization effect (NDVIobs − NDVIadj). Similar to the calculation of the observed ΔNDVImax, 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 NDVIobs and NDVIadj is shown in Supplementary Fig. 7 to highlight temporal nature of this attribution.

### Determining the impact of climate and land use

After the NDVI was scaled to remove the CO2 effect using the relationship from Franks et al.39, a statistical approach was used to attribute the change in NDVIadj to climate (both CV and CC) and LU. We used the recently developed TSS-RESTREND method, which allows vegetation changes due to LU to be separated from those driven by CC and variability19,24,53.

Dryland ecosystems have large natural interannual CV. To separate the impact of climate and LU, previous studies have generally fitted a statistical relationship between climate and vegetation, then used the trends in that relationship or its residuals to quantify LU impacts54,55. However, LU can have both a gradual impact on vegetation through processes such as grazing, which are captured by these methodologies54,56, and abrupt impact through processes such as deforestation, which cause these methods to break down18.

TSS-RESTREND differs from existing dryland trend attribution 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 changes19. To do this TSS-RESTREND incorporates a phenological change detection method57 to identify structural changes in the ecosystem, which manifest as break points in the NDVI time series. We used TSS-RESTREND v2.15, which has been updated to use both precipitation and temperature24, to calculate the Vegetation Climate Relationship (VCR) using the per-pixel optimal precipitation and temperature accumulation periods19,24. TSS-RESTREND was applied to the NDVIadj, with the LU 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 IPCC23, although that approach does not separate CC from CV and also assumes that all dryland plants follow a C3 photosynthetic pathway.

### Separating climate change and climate variability

The TSS-RESTREND method separates the effects of LU from the combined effects of climate (variability and change) using VCR. To separate the effect of CC and CV, the observed climatology (the per-pixel accumulated precipitation and temperature data) was calculated for the period 1962 to 2015. A 20-year leading edge smoothing window was then applied to this observed climatology to remove the interannual CV. The long-term trend caused by CC was determined using the Theil–Sen slope estimator58 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 (NDVICL) and the NDVI due to CV (NDVIcv) were calculated using the observed climatology and detrended climatology, respectively. The difference between NDVICL and NDVICV is the change in NDVImax attributable to CC (NDVICC). The non-parametric Theil–Sen slope estimator and Spearman’s ρ test for significance was then applied to NDVICV and NDVICC to get the change attributable to CV and CC, respectively. The influence of Other Factors (OF), which could not be modeled, was calculated using ({mathrm{{OF}}} = {mathrm{{Obs}}} – ({rm{CO}}_{2} + {rm{LU}} + {rm{CV}} + {rm{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:

$$z_n = frac{{x_n – mu _{mathrm{{obs.}}}}}{{sigma _{mathrm{{obs.}}}}}$$

(3)

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:

$$z_n = frac{{beta times (n – n_0)}}{{sigma _{mathrm{{obs.}}}}}$$

(4)

where β is the per-pixel trend in accumulated precipitation or temperature, n0 is the first year of the analysis (1982). The temperature and precipitation anomaly attributed to CV was calculated using:

$$z_n = frac{{x_{{mathrm{{(adj)}}}n} – mu _{mathrm{{obs.}}}}}{{sigma _{mathrm{{obs.}}}}}$$

(5)

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 Change

In this study, we use the term ACC to refer to both changes in water availability driven by trends in precipitation and increases in temperature14,23, as well as increased water use efficiency (carbon gain per unit of water lost) in response to rising atmospheric CO28,15. The impact of ACC) is calculated using ACC = CO2 + CC. Although we acknowledge that this considering CO2 fertilization 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)

Burrell et al.53 showed that using an ensemble of TSS-RESTREND 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 dryland regions, which can amplify the documented discrepancies between different datasets59,60 We used two 12-member ensembles with matched runs made using TSS-RESTREND 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 CO2, and the second ensemble assumes all plants are C4 with no eCO2 response (see Supplementary Figs. 1 and 2).

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 dryland ecosystems, precipitation controls the amount of foliage cover. As a result, increases in water use efficiency with increasing CO2 will lead to increases in foliage cover in plants that use the C3 photosynthetic pathway8 while plants that use the C4 photosynthetic pathway are not expected to show this response62. Our assumption that C3 plant respond strongly to increased atmospheric CO2, and that C4 plants do not respond, which is consistent with theory and short-term CO2 enrichment experiments62. However, recent surprising results from a long-term CO2 manipulation experiment in a dryland ecosystem have shown much higher responses in C4 plants62. 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 CO2-driven change in NDVI, the Spearman’s rank correlation co-efficient test was used to measure statistical significance for each pixel63. 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 CO2 fertilization effect, the C3 and C4 12-member ensembles included in the Supplementary Material provide an estimate of the upper and lower bounds of CO2 responses (see Supplementary Figs. 1 and 2).

For LU, CV, 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 agreement65. 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 Remapping66 in CDO67. 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 CV, 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 method68.