Global vegetation resilience linked to water availability and variability
Vegetation and land-cover dataTo monitor vegetation at the global scale, we use three datasets: (1) vegetation optical depth (VOD, 0.25°, Ku-Band, daily 1987–201723) (Fig. 1A), (2) AVHRR GIMMSv3g normalized difference vegetation index (NDVI, 1/12°, bi-weekly 1981–201524) (Fig. 1B), and (3) MODIS MOD13 NDVI at 0.05° (16-day, 2000–202125). We correct for spurious values in the NDVI data (e.g., cloud contamination) using the method of Chen et al.43. We resample the VOD data using bi-weekly medians to agree with the NDVI data time sampling.For all three vegetation datasets, we remove seasonality and long-term trends using seasonal trend decomposition by Loess4,44 based on the proposed optimal parameters listed in Cleveland et al.44 (code available on Zenodo45). That is, we use a period of 24 (bi-monthly, 1 year), 47 for the trend smoother (just under 2 years) and 25 for low-pass (just over 1 year). We only use the STL residual—the de-seasoned and de-trended NDVI and VOD time series—in our analysis.To contextualize our understanding of vegetation resilience, we use MODIS MCD12Q1 land cover46 (Fig. 1C) as well as a global average aridity index based on WorldCLIM data31 (Fig. 1D). We exclude from our analysis anthropogenic and non-vegetated landscapes (e.g., permanent snow and ice, desert, urban), as well as any land covers which have changed (e.g., forest to grassland) during the period 2001–2020.Precipitation data and variability metricsTo measure precipitation at the global scale, we rely upon ERA5 data (~30 km, monthly, 1981–2021)33. We process global-scale precipitation metrics using the Google Earth Engine47 platform. We further use the sum of soil moisture from the surface down to 28 cm of depth (first two layers of the ECMWF Integrated Forecasting System soil moisture estimates) to quantify soil moisture means and inter-annual variability33.It is well-documented that vegetation resilience is responsive to the MAP of certain regions1. However, the role of precipitation variability in controlling vegetation resilience has not been well-studied. Here we examine precipitation variability in terms of both intra- and inter-annual patterns. Intra-annual precipitation variability is determined in terms of the Walsh-Lawler Seasonality index32 (Fig. 1D), calculated using monthly data from ERA533.Partly due to the fact that precipitation is non-negative, simple inter-annual variability metrics such as the standard deviation of annual precipitation sums are biased by the absolute precipitation sums; higher precipitation regions have a higher possible range of variability. To limit the influence of MAP, we hence investigate the standard deviation of annual precipitation sums normalized by the MAP, over the period 1981–2021, based on ERA5 data33 (Fig. 1F). We motivate our normalization by MAP with the strong linear relationship between MAP and MAP standard deviation (Supplementary Fig. S2). We further confirm our discovered relationships (Fig. 5) using only those regions where MAP was between the 40 and 60th percentile of MAP for a given land cover (Supplementary Figs. S11,S12). This serves as an additional check that our normalization of MAP standard deviation by MAP does not bias the inferred relationship between vegetation resilience and precipitation variability. Similarly, we generate a normalized inter-annual soil moisture variability by normalizing year-on-year soil moisture standard deviation (Supplementary Fig. S8) by long-term mean soil moisture (Supplementary Fig. S5).Empirical resilience estimationResilience is defined as the ability of a system to recover from perturbations, and can be quantified empirically by the speed of recovery to the previous state16,17. To measure resilience on the global scale, we employ a recently introduced methodology4 which we will briefly summarize in the following.We first identify sharp transitions in the vegetation time series using an 18-point (9 month) moving window to define local slopes throughout the time series48. We then identify slopes above the 99th percentile, and define connected regions as individual perturbations. The highest peak (largest instantaneous slope) within each connected region is then labeled as an individual disturbance.The employed approach does not delineate every rapid transition in a time series due to our reliance on percentiles; our dataset will be inherently biased towards the largest transitions. Furthermore, the same transitions are not guaranteed to be captured for both NDVI and VOD data in each location, as the percentiles will naturally vary between the datasets. Finally, our method will in some cases produce false positives, especially in cases where a given time series does not have any significant rapid transitions. To limit the influence of false positives on our results, we discard any perturbations where the time series does not drop significantly, and where the period before and after a given transition does not pass a two-sample Kolmogorov–Smirnov test4.Finally, using our global set of time-series transitions, we can identify each local vegetation (NDVI or VOD) minima, and use the five following years of data to fit an exponential function to the residual time series, assuming that the recovery after a perturbation to a vegetation state x0 follows approximately the equation$$x(t),approx ,{x}_{0}{e}^{rt}$$
(1)
where x(t) denotes the vegetation state at time t after the perturbation. Negative r indicates that the vegetation system will return to the original stable state at rate ∣r∣. For positive r, the initial perturbation would be amplified, suggesting a non-resilient vegetation state. Our empirical recovery rates are defined as the fitted exponent r, obtained for each detected transition in the NDVI and VOD residual time series. We finally use the coefficient of determination R2 to remove instances where the fitted exponential poorly matches the underlying data4.For the empirical estimate of the restoring rate obtained from fitting an exponential to the recovery after an abrupt negative deviation of VOD or NDVI, abrupt changes in the mean state induced by changing sensors rather than an actual vegetation shift may impact the results. However, all datasets used here are tightly cross-calibrated to eliminate mean-shifts when new instruments are introduced23,24. It is therefore unlikely that changes in the instrumentation of the various datasets unduly influence our empirical estimates of λ.Dynamical system metrics of resilienceThe lag-one autocorrelation (AC1) has previously been proposed to measure the stability of real-world dynamical systems in general, and the resilience of vegetation systems in particular1,19,20,21,49. Based on the concept of critical slowing down, the AC1 has, together with the variance, also been suggested as an early-warning indicator for forthcoming critical transitions50,51. Mathematically, the suitability of the variance and AC1 as resilience measures and early-warning indicators can be motivated as follows4,52,53. First, linearize the system around a given stable state x*:$$dbar{x}=lambda bar{x}dt+sigma dW$$
(2)
for (bar{x}: !!=x-{x}^{*}), assuming a Wiener Process W with standard deviation σ. The dynamics are stable for λ More
