Global vegetation resilience linked to water availability and variability

Vegetation and land-cover data

To monitor vegetation at the global scale, we use three datasets: (1) vegetation optical depth (VOD, 0.25°, Ku-Band, daily 1987–2017²³) (Fig. 1A), (2) AVHRR GIMMSv3g normalized difference vegetation index (NDVI, 1/12°, bi-weekly 1981–2015²⁴) (Fig. 1B), and (3) MODIS MOD13 NDVI at 0.05° (16-day, 2000–2021²⁵). We correct for spurious values in the NDVI data (e.g., cloud contamination) using the method of Chen et al.⁴³. 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 Loess^4,44 based on the proposed optimal parameters listed in Cleveland et al.⁴⁴ (code available on Zenodo⁴⁵). 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 cover⁴⁶ (Fig. 1C) as well as a global average aridity index based on WorldCLIM data³¹ (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 metrics

To measure precipitation at the global scale, we rely upon ERA5 data (~30 km, monthly, 1981–2021)³³. We process global-scale precipitation metrics using the Google Earth Engine⁴⁷ 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 variability³³.

It is well-documented that vegetation resilience is responsive to the MAP of certain regions¹. 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 index³² (Fig. 1D), calculated using monthly data from ERA5³³.

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 data³³ (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 estimation

Resilience 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 state^16,17. To measure resilience on the global scale, we employ a recently introduced methodology⁴ 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 series⁴⁸. 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 test⁴.

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 x₀ follows approximately the equation

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 R² to remove instances where the fitted exponential poorly matches the underlying data⁴.

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 introduced^23,24. It is therefore unlikely that changes in the instrumentation of the various datasets unduly influence our empirical estimates of λ.

Dynamical system metrics of resilience

The 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 particular^{1,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 transitions^50,51. Mathematically, the suitability of the variance and AC1 as resilience measures and early-warning indicators can be motivated as follows^4,52,53. First, linearize the system around a given stable state x^*:

for (bar{x}: !!=x-{x}^{*}), assuming a Wiener Process W with standard deviation σ. The dynamics are stable for λ < 0 and unstable otherwise. Upon discretizing the resulting Ornstein-Uhlenbeck process into time steps of width Δt, the variance and AC1 of the resulting order-one autoregressive process are then related to the restoring rate λ via⁵⁴:

for the variance, and

for the AC1. Hence, the closer λ is to zero, the larger the AC1 and variance, corresponding to lower stability. Note that theory⁵⁵ suggests that the recovery rate r is equal to the restoring rate λ. An empirical global confirmation for this relationship for has recently been demonstrated based on both NDVI and VOD data⁴. We compare the dependence of the recovery rate r and two different theoretical estimates of the restoring rate λ—obtained via inverting the above equations for the variance and AC1—with respect to their dependence on water availability and its variability.

It is important to note that combining data from different sensors with varying signal-to-noise ratios (e.g., VOD, AVHRR NDVI) can bias estimates of temporal changes in resilience indicators because the higher-order statistics of the resulting time series are not homogeneous^4,26. In the present work, however, we do not investigate temporal trends via estimating resilience indicators in sliding windows (as in refs. ^4,26), but rather estimate resilience indicators for the full available time series. This excludes the possibility of systematic biases in our AC1- and variance-based estimates of the restoring rate λ. In principle, the combination of different sensors might lead to larger uncertainties in the estimates of λ: combining different sensor data leads to temporally varying (yet spatially homogeneous) effects on the AC1 and variance and may therefore lead to a wider spread—but not a bias—in the resulting estimates of λ.

Binning and significance testing

The direction and magnitude of our discovered relationships (e.g., Figs. 3–5, Supplementary Figs. S4–S12) will be to some degree controlled by the choice of bin sizes. We tested three bin sizes for each different variable (Supplementary Figs. S13,S14). We also imposed the conditions that there were at least 50 measurements in each bin to form a proper median, and that we only report those relationships which cover ten or more bins.

To better constrain the relationship between each driving variable and resilience, we use the non-parametric Kendall-Tau test⁵⁶. Kendall-Tau statistics are calculated over each set of binned medians, as well as using a Monte-Carlo approach. Over 1000 iterations, we choose one random point from each bin and recalculate the Kendall-Tau statistics. These 1000 surrogates are displayed as box plots in Figs. 3–5. Note that the Kendall-Tau value of the median line will almost always be larger than the median of the 1000 Kendall-Tau values resulting from the surrogates due to the smoothing inherent in taking binned medians. That is, the binned medians represent a smooth and almost monotonic line with fewer jumps, while the 1000 surrogates will have strong fluctuations from one bin to the next, leading to overall lower Kendall-Tau values. The fraction of Kendall-Tau statistics which share a sign with the Kendall-Tau of the median line is also reported on Supplementary Figs. S4–S14.

We find that our Kendall-Tau statistics are robust against changing bin sizes (Supplementary Figs. S13,S14), where the direction of trends does not change from those reported in Fig. 3. The magnitude of the Kendall-Tau statistic—as well as p-values—shift with different bin sizes, with smaller bin sizes typically resulting in more robust trends. Changes in bin size do not have a strong impact upon our data interpretations or conclusions.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com