Satellite-Observed Global Terrestrial Vegetation Production in Response to Water Availability

Abstract

Water stress is one of the primary environmental factors that limits terrestrial ecosystems’ productivity. Hense, the way to quantify gobal vegetation productivity’s vulnerability under water stress and reveal its seasonal dynamics in response to drought is of great significance in mitigating and adapting to global changes. Here, we estimated monthly gross primary productivity (GPP) first based on light-use efficiency (LUE) models for 1982–2015. GPP’s response time to water availability can be determined by correlating the monthly GPP series with the multiple timescale Standardized Precipitation Evapotranspiration Index (SPEI). Thereafter, we developed an optimal bivariate probabilistic model to derive the vegetation productivity loss probabilities under different drought scenarios using the copula method. The results showed that LUE models have a good fit and estimate GPP well (R² exceeded 0.7). GPP is expected to decrease in 71.91% of the global land vegetation area because of increases in radiation and temperature and decreases in soil moisture during drought periods. Largely, we found that vegetation productivity and water availability are correlated positively globally. The vegetation productivity in arid and semiarid areas depends considerably upon water availability compared to that in humid and semi-humid areas. Weak drought resistance often characterizes the land cover types that water availability influences more. In addition, under the scenario of the same level of GPP damage with different drought degrees, as droughts increase in severity, GPP loss probabilities increase as well. Further, under the same drought severity with different levels of GPP damage, drought’s effect on GPP loss probabilities weaken gradually as the GPP damage level increaes. Similar patterns were observed in different seasons. Our results showed that arid and semiarid areas have higher conditional probabilities of vegetation productivity losses under different drought scenarios.

1. Introduction

The Intergovernmental Panel on Climate Change (IPCC) report shows that the world will continue to warm in the 21st century [1]. This indicates that as the future temperature rises, the occurrence and frequency of extreme weather and climate, such as heatwaves and droughts, will increase rapidly [2]. This trend and its associated adverse effects on natural and social ecosystems are expected to increase further [3]. Terrestrial ecosystems’ gross primary productivity (GPP) is the largest component of global terrestrial carbon flux, and slight fluctuations in GPP have a significant influence on atmospheric CO₂ concentration [4]. On a global scale, particularly in arid and semiarid regions, water stress is the primary environmental factor that limits terrestrial ecosystems’ productivity [5]. Hence, studying drought’s effects on their productivity has become an important priority in global change research.

The terrestrial ecosystem’s productivity characterizes its quality, it is the key to controlling the ecosystem’s carbon, water, energy and nutrient cycles, and plays a decisive role in regulating its function and controlling its carbon source and sink [6]. At the same time, ecosystem productivity is also the basis of a variety of ecosystem services [7] and is the source of food, fiber, wood, livestock grains and biofuels in human society [8]. As the main outcome and important manifestation of global climate change, drought not only affects photosynthesis and ecosystem productivity directly [9], but it can also affect terrestrial ecosystem productivity in directly by changing the intensity and frequency of other forms of disturbance, such as increasing the frequency and intensity of fires, plant mortality and the occurrence of pests and diseases [10]. In recent years, many studies have investigated drought’s effects on terrestrial ecosystems’ productivity at different spatio-temporal scales based on long-term observations from ecological research sites, field experiments, eddy covariance measurements, large-scale satellite measurements and model simulations [11,12,13]. Satellites data used commonly to evaluate vegetation productivity include primarily the Moderate Resolution Imaging Spectroradiometer-Fraction of Photosynthetically Active Radiation (MODIS-FPAR) and the Global Inventory Modeling and Mapping Studies-Fraction of Photosynthetically Active Radiation (AVHRR GIMMS-FPAR), remotely-sensed data for atmosphere state, MODIS Enhanced Vegetation Index, leaf area index, spatially distributed CO₂ concentrations, and canopy information. For example, Reichstein et al. [11] found that the drought stress in Europe in summer 2003 caused the reduction of ecosystem productivity by jointing flux tower, remote sensing (including MODSI and AVHRR-GIMMS FPAR) and modelling datasets. Zhang et al. [14] emphasized that the drought in Southwestern China in spring 2010 reduced primary productivity by using primary productivity products derived from MODIS. Stocker et al. [15] used four remote sensing data-driven models and demenstrated that drought’s effects on terrestrial primary production were underestimated. Previous studies have shown that different vegetation types and environmental conditions (such as climate and soil factors) largely determine biological communities’ resistance and ability to recover from drought stress [16].

Although there is increasing recognition of drought’s adverse effects on terrestrial vegetation productivity, few studies have focused on the seasonal dynamics in global productivity in response to drought across various climate zones and land biomes. Moreover, studies of drought’s lag effect on terrestrial vegetation productivity and vegetation’s responses to drought at various time-scales on a global scale are scare. The Standardized Precipitation Evapotranspiration Index (SPEI) quantifies different drought types based on the water balance, and is used widely on global and regional scales because of its multiple time scales and spatial comparability [17]. Therefore, determining the corresponding time scale when GPP and SPEI have the highest correlation based upon different time scales is of great significance in understanding the variations in different terrestrial ecosystems’ productivity depending uopn drought resistance. In addition, we give greater attention here to the way to quantify the probability of varying degrees of damage to vegetation productivity under predictable drought scenarios. Copulas can be used to couple multiple variables and construct their joint probability distribution by considering their correlations [18], which can derive the vegetation productivity loss distributions conditioned on any drought scenario. For example, Madadgar et al. [19] estimated the probability of drought in agricultural production in Australia based on a copula model, while Fang et al. [20] constructed a bivariate probabilistic framework to assess vegetation vulnerability and map drought-prone ecosystems in China’s Loess Plateau.

Given ecosystem productivity’s key role in affecting the global carbon budget and the supply of ecosystem services, in the context of global warming, the rapid increase in atmospheric carbon dioxide concentration, and increasing human food production and resource demand, estimating ecosystem productivity accurately is important to assess drought’s effects on terrestrial ecosystems precisely [21]. Three main methods are used to estimate ecosystems’ productivity at the regional and global scales—remote sensing-based light-use efficiency (LUE) models [22], machine learning algorithms [23], and process-based physiological and ecological models [24]. Remote sensing observations can provide information continuously on land surface features that affect ecosystem productivity, such as ecosystem structure, vegetation phenology, biomass, soil moisture, and land cover over a large area. The LUE model relies on remote sensing observation data to estimate ecosystem productivity through the utilization rate of photosynthetically active radiation in the ecosystem [25]. We gave special attention to determining water stress and produced a data set of monthly LUE GPPs at the global scale based on various representations of water stress: Water stress based on vapor-pressure deficit, humidity deficit and root-zone soil-water content separately.

Our study had two primary objectives. The first was to determine the seasonal dynamics of global vegetation productivity’s response to drought across various climate zones and land biomes, and the second was to quantify the probability of varying degrees of damage to vegetation productivity under predictable drought scenarios. First, we used a LUE model to produce a dataset of global monthly satellite-observed GPPs, which are designed primarily to determinate water stress. Second, we explored vegetation productivity’s spatio-temporal evolution of dependence and response time to water availability by correlating monthly GPP series with the multiple timescale SPEIs (1- to 24-months). Finally, we develpoed an optimal bivariate probabilistic model in each pixel of the globe using the copula method, which was used to derive the vegetation productivity loss probabilities under different drought scenarios.

2. Materials and Methods

2.1. Data

2.1.1. Remote Sensing Data

The actual evapotranspiration, potential evapotranspiration and soil moisture (including surface soil and root-zone soil moisture) data were derived from the Global Land Evaporation Amsterdam Model (GLEAM), which is designed to be determined by remote sensing observations only [26]. The GLEAM_V3.3a dataset has a 0.25° spatial resolution at a monthly time step that spans the 39-year period from 1980–2018 (www.gleam.eu (accessed on 25 March 2021)). Fraction of Photosynthetically Active Radiation (FPAR)3g data were generated from AVHRR GIMMS NDVI3g using an Artificial Neural Network (ANN) derived model, with a spatial resolution of 1/12° and a 15 days time interval [27]. We used the MODIS MCD12C1 product, which contains yearly, worldwide distributions of 17 land cover types from 2001–2017 at a spatial resolution of 0.05° (https://lpdaac.usgs.gov/ (accessed on 10 December 2020)) [28]. In this study, barren land, water bodies, permanent snow, and ice were removed.

2.1.2. Meteorological data and FLUXNET data

Downward shortwave radiation, air temperature, specific humidity and air pressure were obtained from the CRU-NCEPV6.1 dataset with a spatial resolution of 0.5° and a 6-h time interval. We used gross primary productivity (GPP) data from the FLUXNET2015 dataset to compare flux site measurement of GPP and LUE GPP (https://fluxnet.org/data/fluxnet2015-dataset/ (accessed on 1 December 2020)).

2.2. Drought Index

SPEI can be calculated by fitting the difference between precipitation and potential evapotranspiration based on a log-logistic probability distribution, which combines the multi-temporal characteristic of standardized precipitation index and the sensitivity of palmer drought severity index to changes in evaporation demand. A global gridded dataset of the SPEI was calculated using monthly precipitation and potential evapotranspiration from the CRU TS3.24.01 dataset [29]. This dataset is generated at a 0.5° spatial resolution and at a monthly time step and covers the period from 1902 to 2015. We used the SPEI dataset from 1982 to 2015 in this study, which contains timescales from 1- to 24-months. A 6-month SPEI value is formulated by the cumulative water deficit or surplus from five months before to the current month. Three drought scenarios including moderate (−1.5 < SPEI ≤ 1), severe (−2 < SPEI ≤ −1.5), and extreme (SPEI ≤ −2) can be differentiated using a set of SPEI thresholds [29]. The yearly aridity index (AI), which has a spatial resolution of 0.5°, was calculated with Feng and Fu’s formula [30]. (Response: we cited the citation 31 in Section 2.3) According to this dataset, global land can be classified into arid (0.05 ≤ AI < 0.2), semiarid (0.2 ≤ AI < 0.5), semihumid (0.5 ≤ AI < 0.65) and humid zones (AI ≥ 0.65).

2.3. LUE GPP

The global GPP for the 1982–2015 period was calculated based on the LUE equation, as follows [31]:

$$GPP={\epsilon }_{max}\times SOL\times 0.45\times fPAR\times f\left(T\right)\times f\left(W\right),$$

(1)

in which ${\epsilon }_{max}$ is vegetation type related maximum light use efficiency, $SOL$ is downward solar shortwave radiation, and $fPAR$ represents the fraction of photosynthetically active radiation vegetation absorbs. $ft$ and $fw$ are temperature and water stress respectively to plant productivity. We adopted the meteorological data from the CRU-NCEPV6.1 dataset, which includes downward shortwave radiation, air temperature, specific humidity and air pressure, and so forth. We used $fPAR$ Zhu et al. [27] developed. We calculated the temperature stress with Equation (2), as follows:

$$ft=\frac{\left(T-T_{MIN}\right)\left(T-T_{MAX}\right)}{\left[\left(T-T_{MIN}\right)\left(T-T_{MAX}\right)\right]-{\left(T-T_{opt}\right)}^2},$$

(2)

in which $T_{MIN}$ and $T_{MAX}$ are the minimum and maximum temperature for plant’s photosynthesis, with the specific parameters can be found in Zhang et al. [32].

In this study, we focused particularly on determining water stress ($fw$) in the LUE equation. The MODIS-GPP algorithm [33] calculates $fw$ with the daily vapor pressure deficit (VPD) on consideration that leaf stomatal conductance is mainly restricted by cold stress and VPD, named VPDMOD expressed as in Equation (3):

$$fw=f\left(VPD\right)=\{\begin{array}{c}1 VPD\le VP{D}_{open}\\ \frac{VP{D}_{close}-VPD}{VP{D}_{close}-VP{D}_{open}} VP{D}_{open}<VPD<VP{D}_{close}\\ 0 VPD\ge VP{D}_{close}\end{array},$$

(3)

in which $VP{D}_{close}$ represents VPD that induces full stomatal closure, while $VP{D}_{open}$ represents VPD that results in full stomatal opening, all parameters are set based on the vegetation type.

The GLO-PEM algorithm [34] considers each vegetation type’s optimum growth temperature ($T_{opt}$). The algorithm calculates water stress with vegetation type, independently with Equation (4) named VPDGLO:

$$fw=f\left(VPD\right)=\{\begin{array}{c}1-0.05{\delta }_q 0<{\delta }_q\le 15\\ 0.25 {\delta }_q>15\end{array},$$

(4)

in which ${\delta }_q$ is the specific humidity deficit (g/kg), the difference between saturated and actual specific humidity.

The CASA model provides a framework that considers water stress with a multidisciplinary algorithm [35], shown in Equation (5), as follows:

$$\begin{gathered} fw=f\left(VPD\right)\times W_s \\ {W}_s=\{\begin{array}{c}1 VPD only model\\ 0.5+ETR VPD+ETR model\\ 0.5 + K_S VPD+SM model\end{array}, \end{gathered}$$

(5)

in which $ETR$ is the ratio of actual evapotranspiration to potential evapotranspiration. According to FAO 56 [36], $ETR$ can be equal to $K_S$, which is calculated based on root zone soil moisture (hereinafter $SM$) and soil and vegetation properties, written as Equations (6) and (7):

$$K_S=\frac{SM-WP}{\left(1-p\right)\times TAW}$$

(6)

$$p=p_{std}+0.04\times \left(5-E{T}_0\right),$$

(7)

in which WP and $TAW$ represent soil wilting point and total available soil water derived from the International Geosphere-Biosphere Programme Data and Information System (IGMP-DIS) dataset [37] and the World Inventory of Soil property Estimates (WISE) derived soil properties [38]. $E{T}_0$ is calculated with Hagreaves equation [39] while $p_{std}$ is a vegetation type-specific depletion factor (allowable extent of soil water deficit before water stress) at $E{T}_0$ of 5 mm per day. Actual evapotranspiration, potential evapotranspiration, and root-zone soil moisture data were all acquired from the GLEAM dataset.

In this study, we provided a dataset of monthly plant productivity that spans the 34-year period from 1982–2015 at the global scale, composed by five LUE GPPs derived from the multidisciplinary framework in Equation (5). These five GPPs represent different considerations of water stress in the LUE equation, the VPD only model, VPD_GLO-ETR model, VPD_MOD-ETR model, VPD_GLO-SM model, and VPD_MOD-SM model. Because the LUE equation estimate GPP from multiple satellite data, for example, NDVI, FPAR, and GLEAM soil water, it is considered satellite-observed GPP.

2.4. Determining LUE GPP’s Response Time to Water Availability

To screen out vegetation productivity’s response time to water availability, a monthly GPP series from 1982 to 2015 was correlated with the multiple timescale SPEI (1- to 24-months). For the i-th month in a year, Spearman’s correlation analysis was performed between the monthly GPP and SPEI series at different timescales as follows

$$R_j^i=corr\left(GP{P}^i,SPE{I}_j^i\right) i=1, 2 ,\dots ,12, j=1, 2,\dots ,24,$$

(8)

in which R is the correlation coefficient and a timescale j denotes the cumulative water deficit or surplus from j-1 months before to the current month. Finally, a timescale able to maximize the GPP-SPEI association, is retained as the vegetation productivity response time (VPRT) to water variations for the i-th month.

$$VPR{T}^i=max\left\{R_j^i\right\} j=1, 2, \dots , 24.$$

(9)

2.5. Copulas

We chose the extreme value (EV), generalized extreme value (GEV), exponential (EXP), gamma (GAM), Poisson (POISS), normal (NOR), generalized Pareto (GP) and Weibull (WBL) distribution functions to fit the monthly GPP series (Table S1). The SPEI series is a statistic distributed normally and accordingly, the normal distribution is employed as its marginal. The probability distribution functions of monthly GPP and SPEI series are defined as F_GPP (gpp) and F_SPEI (spei), respectively.

We constructed a joint distribution function in each pixel of the globe in every month based on monthly GPP and SPEI series. Considering F_GPP (gpp) and F_SPEI (spei), these two marginal distribution’s joint distribution function can be defined as a copula C:

$$F_{GPP,SPEI}\left(gpp,spei\right)=C\left(F_{GPP}\left(gpp\right),F_{SPEI}\left(spei\right)\right).$$

(10)

In this study, we used the methods of maximum likelihood and inference functions for margins to estimate the fitted parameters in F_GPP (gpp), F_SPEI (spei) and F_GPP,_SPEI (gpp, spei) [40]. The bivariate distribution functions used commonly, including the Frank, Clayton, Gumbel, t and Normal copulas, were constructed to couple FGPP (gpp)) and FSPEI (spei) in each pixel of the globe (as shown in Table S2).

Next, we obtained the optimal joint distribution function of F_GPP (gpp) and F_SPEI (spei) as follows by referring to the method of the goodness-of-fit test Zhang et al. [41] provided.

$$F_{GPP,SPEI}\left(gpp,spei\right)=P\left(GPP\le gpp,SPEI\le spei\right)={\int }_{-\infty }^{gpp}{\int }_{-\infty }^{spei}f\left(gpp,spei\right)d_u{d}_v.$$

(11)

Ultimately, we generated a set of conditional probabilities of vegetation productivity losses for diverse drought conditions with a bivariate statistical framework. Larger values of vegetation productivity loss probabilities under the same scenario imply greater ecosystem vulnerability, and such ecosystems are categorized thereby as the drought-prone class. Given multiple drought scenarios, conditional probabilities of vegetation productivity lower than different percentiles (e.g., the 10th, 20th, 30th, 40th percentiles were considered in the study) are derived using copula-based joint distribution and the conditional distribution formulas. The conditional probabilities of different percentiles of vegetation productivity loss under moderate, severe and extreme drought scenarios can be calculated by the following formula.

$$P\left(GPP<gpp| -1.5<SPEI\le -1\right)=\frac{F_{GPP, SPEI}\left(gpp,-1\right)-F_{GPP, SPEI}\left(gpp, -1.5\right)}{F_{SPEI}\left(-1\right)-F_{SPEI}\left(-1.5\right)}$$

(12)

$$P\left(GPP<gpp|-2<SPEI\le -1.5\right)=\frac{F_{GPP, SPEI}\left(gpp,-1.5\right)-F_{GPP, SPEI}\left(gpp, -2\right)}{F_{SPEI}\left(-1.5\right)-F_{SPEI}\left(-2\right)}$$

(13)

$$P\left(GPP<gpp| SPEI\le -2\right)=\frac{F_{GPP, SPEI}\left(gpp,-2\right)}{F_{SPEI}\left(-2\right)}.$$

(14)

2.6. Statistical Tests

We used the linear regression method to calculate the GPP’s mean change rate over the entire temporal domain. A two-tailed t-test and the Mann-Kendall test were used to examine the significance of the trend, while the Kolmogorov–Smirnov test was used to establish the optimal marginal distribution function of GPP (p < 0.05). In addition, we adopted several goodness-of-fit measures to evaluate different copula models’ performance, including the squared Euclidean distance (SED) between the theoretical copula and the empirical copula, the root mean squared error (RMSE) and the Akaike information criterion (AIC).

3. Results

3.1. Validation of LUE GPPs’ Accuracy, Dynamics Trends and Drought’s Effect on GPP

In this study, special attention was given to determining water stress (including atmosphere vapor pressure deficit, soil moisture content and the humidity deficit) in the LUE model to estimate the global GPP. By comparing with the FLUXNET GPP site data, we found that different LUE models have a good fit in estimation of LUE GPP. The fitted R² of VPD_GLO-SM, VPD_MOD-SM, VPD_GLO-ETR, and VPD_MOD-ETR were 0.7739, 0.7399, 0.7427, 0.7459 and 0.7628, respectively (

Table 1. Regression statistics for different light-use efficiency models between each land-use type’s measured FLUXNET GPP and LUE GPPs.

		Light-Use Efficiency Models
Type	Statistics	VPD only	VPDGLO-SM	VPDMOD-SM	VPDGLO-ETR	VPDMOD-ETR
DBF (N: 1693)	R2	0.7822	0.7467	0.7610	0.7275	0.7454
	RMSE	2.1887	2.3839	2.3068	2.4894	2.3932
EBF (N: 857)	R2	0.6388	0.3798	0.3663	0.5599	0.6662
	RMSE	1.9865	2.8458	2.8692	2.5501	1.9500
ENF (N: 2968)	R2	0.7706	0.7833	0.7733	0.7671	0.7590
	RMSE	1.6990	1.6399	1.6919	1.7278	1.7697
MF (N: 863)	R2	0.8480	0.8479	0.8484	0.8486	0.8496
	RMSE	1.3900	1.3825	1.3877	1.3873	1.3904
WET (N: 541)	R2	0.6860	0.5871	0.6297	0.5953	0.6406
	RMSE	2.0777	2.5013	2.3196	2.4670	2.2742
CSH/OSH (N: 317)	R2	0.6563	0.4718	0.4369	0.6829	0.6609
	RMSE	1.3422	1.7781	1.8739	1.3185	1.3611
WSA (N: 545)	R2	0.8103	0.7806	0.7947	0.7784	0.8001
	RMSE	1.1185	1.3406	1.2351	1.3736	1.2354
SAV (N: 427)	R2	0.6776	0.6511	0.6517	0.6749	0.6746
	RMSE	1.3769	1.6252	1.6233	1.4901	1.4919
GRA (N: 1767)	R2	0.7667	0.7727	0.7730	0.7754	0.7763
	RMSE	1.9214	1.9175	1.9088	1.9137	1.9012
CRO (N: 1392)	R2	0.5290	0.4686	0.5047	0.4564	0.4953
	RMSE	3.5701	/
All sites except cropland (N: 9978)	R2	0.7739	0.7399	0.7427	0.7459	0.7628
	RMSE	1.8089	1.9817	1.9679	1.9739	1.8855

Note:The sites are deciduous broadleaf forest (DBF), evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), mixed forest (MF), wetlands (WET), closed shrublands/open shrublands (CSH/OSH), woody savanna (WSA), savanna (SAV), grassland (GRA) and cropland (CRO). VPD only, VPD_GLO-SM, VPD_MOD-SM, VPD_GLO-ETR and VPD_MOD-ETR stand for the moisture stress algorithms. N represents the total number of data points; R² is the R²e values of linear regressions with intercept. RMSE represents root mean square error. For each land-use type, the corresponding “best algorithm” is highlighted in bold.

). The average goodness of fit of multiple models reached 0.78 (Figure S2). Overall, LUE GPP had a greater response to the atmospheric vapor pressure deficit, followed by soil water content, and the humidity deficit. The LUE GPP showed greater sensitivity to soil moisture content than the vapor pressure deficit only in evergreen needleleaf forest (ENF), mixed forest (MF) (R² = 0.7833 and 0.8479, respectively, RMSE = 1.6399 and 1.3825, respectively). The annual average spatial distribution of GPP in multiple models is shown in Figure S1. From 1982 to 2015, the global average annual GPP of terrestrial vegetation continued to increase at a mean rate of 0.134 Pg C a ⁻¹ (p < 0.001), but its growth rate declined after the mid-1990s (Figure 1a). From the spatial distribution of the mean annual GPP change trend, it can be seen that the global mean annual GPP still had a largely increasing trend, and the significant increase area accounted for 36.83% of the global terrestrial vegetation area, located primarily in the mid-high latitudes of the northern hemisphere, India, Southeastern China, Southeastern South America, Northern and Southern Africa, and Eastern Australia. 11.62% of the regional GPP showed a significant downward trend, largely in Eastern Brazil and Central and Southern Africa (Figure 1b).

Taking SPEI-12 as an example, a drought episode with an SPEI value less than −1 and no less than 3 consecutive months was considered a drought event. We counted the total frequency of droughts from 1982 to 2015 in each grid around the world. It can be seen that the global drought-prone areas are concentrated primarily in Central Russia, Southeast Asia, the Mediterranean coast and most parts of Africa, and the frequency of droughts has increased more than 10 times (Figure 2a). We divided the entire time period from 1982 to 2015 into a dry period and a non-dry period, and compared the changes in LUE GPP during the dry period relative to the non-dry period. It can be seen that in the high latitudes of the Northern hemisphere, Eastern Asia, the Amazon basin and Central Africa, vegetation productivity under drought conditions increased (accounting for 28.09% of the total global land vegetation area), while in 71.91% of the total land vegetation area in the world, the occurrence of drought led to a decrease in GPP (Figure 2b). To reveal the cause of this phenomenon further, we considered the effect of climate factors (including temperature, solar radiation and soil moisture) on vegetation productivity. As shown in Figure 3, temperature and solar radiation play important roles in the changes in GPP in the high latitudes of the Northern hemisphere, Eastern Asia, the Amazon basin and Central Africa during drought periods, while in other regions, soil moisture limits GPP primarily. In arid and semiarid areas, GPP is correlated negatively with temperature and solar radiation, while in humid and semi-humid areas, the converse is true. GPP is correlated positively with soil moisture globally.

3.2. Spatio-Temporal Dynamics of Vegetation Productivity’s Dependence on Water Availability

The monthly LUE GPP series was correlated with the SPEI of multiple time scales (1–24 months) to determine the extent to which water resources determine global vegetation productivity changes. As shown in Figure 4a, vegetation productivity and water availability are largely correlated positively globally, indicating that vegetation productivity tends to increase (decrease) as water supply increases (decreases). In 47.30% of the global terrestrial ecosystem, the correlation coefficient between GPP and SPEI can exceed 0.6. Particularly in the Southwest United States, Southeastern Argentina, South Africa, Central Asia, and Australia, the correlation coefficient between GPP and SPEI can exceed 0.8, indicating that vegetation productivity in these areas depends more upon water supply. Seasonal changes also affect the correlation between GPP and SPEI. As Figure S3a–d shows, the regions with a high correlation between GPP and SPEI in the United States and Central Asia were located primarily in the Southwest of these two regions in March–May, and from June to August, the scope and intensity of GPP that SPEI affected in these two regions reached the maximum (the correlation coefficient was 0.6 or higher). As the rainy season in the Northern hemisphere fades, the correlation between GPP and SPEI in the United States and Central Asia weakens gradually, and from December to February of the following year, the correlation coefficient between GPP and SPEI in this region fell below 0.4. Vegetation productivity and water availability’s dependence in the high latitudes of the Northern hemisphere is relatively stable, and remains largely between 0.2 and 0.4, while the correlation between GPP and SPEI in Australia is strong throughout the year.

As shown in Figure 5a, an analysis of the changes in the correlation between GPP and SPEI in various climate zones indicated that as the climatic conditions become gradually humid, the dependence between vegetation productivity and water availability continues to decrease (the correlation coefficient between GPP and SPEI declined from 0.76 to 0.47). Different seasons showed consistent characteristics, indicating that the vegetation productivity in arid and semiarid areas is more sensitive to changes in water availability than in humid and semi-humid areas. By comparing the maximum correlation coefficients of GPP and SPEI under different land cover types (Figure 6a), it can be seen that the productivity of GRA, SAV and DBF has a greater correlation with water availability (the correlation coefficients are 0.69, 0.64 and 0.61 respectively), followed by CRO, MF, DNF, ENF, OS and EBF (correlation coefficients of 0.57, 0.55, 0.53, 0.51, 0.50, and 0.50 in order). WSA and WET’s productivity had the weakest correlation with water resources. Our results showed that various land cover types have different adaptation strategies to the increase and loss of water resources.

3.3. Spatio-Temporal Dynamics of Vegetation Productivity’s Response Time to Water Availability

We analyzed the temporal and spatial dynamics of vegetation productivity’s response time of to water availability changes further. Productivity’s response time to water availability is defined as timescales over which the maximum GPP-SPEI correlations are recorded. It can be seen from the spatial distribution of productivity’s response time to water availability (Figure 4b) that the response time of 56.8% of the global terrestrial ecosystems to water resources is based primarily on short-term and medium-term time scales. In Canada, Northwestern Colombia, Eastern Russia, and Northeastern and Southeastern China, the response time of vegetation productivity to SPEI is relatively long, typically more than 12 months. The longer the response time to water availability changes, the stronger the ability of vegetation productivity in these areas to withstand long-term water shortages. Seasonal changes also affect the length of time in which GPP responds to SPEI in different regions. It can be seen from Figure S3e–h that 34.92% of the global terrestrial ecosystem GPP requires a long time to respond to SPEI during the period from March to May and is concentrated largely in the mid-high latitudes of the Northern hemisphere. With the passage of seasons, the proportion of GPPs in these high latitude regions with a longer response time to SPEI decrease gradually from June to November. From December to November of the following year in Australia, GPP’s response time to SPEI changes from a short-term time scale (less than 6 months) to a medium- and long-term time scale (6–9 months).

As shown in Figure 5b, an analysis of the changes GPP’s in the response time to SPEI in various climate zones indicates that as the climatic conditions become gradually humid, the response time of vegetation productivity to changes in water availability increased (GPP’s mean response time to SPEI extended from 3.9 months to 8.9 months). The characteristics were consistent in different seasons, indicating that vegetation’s productivity capacity to withstand long-term water shortages is weaker in arid and semiarid areas than in humid and semi-humid areas. By comparing GPP’s response time to SPEI under different land cover types (Figure 6b), it can be seen that DNF has the strongest ability to resist long-term water shortages (the mean response time reached 13.6 months), followed by OS, EBF, WET, CRO, ENF, WSA and MF (mean response time of 9.9, 8.5, 8.4, 8.1, 7.9, 7.8 and 7.5 months, respectively). SAV, GRA, and DBF have poor ability to resist and mitigate drought (mean response time of 6.1, 6, and 4.9 months, respectively). It is worth noting that the land cover types that are more relevant to water availability are often accompanied by weak drought resistance.

3.4. Vegetation Productivity Loss Probability under Different Drought Scenarios

We first fitted GPP’s optimal marginal distribution during different months based on the Kolmogorov–Smirnov test. According to the criterion that the smaller the SED, RMSE, and AIC values, the better the copula function’s fit effect, we fitted the optimal copula function with the GPP and SPEI in each grid. We set GPP’s damage degree to four levels, the 10th, 20th, 30th, and 40th percentile values of GPP in different seasons. Figure 7 shows the spatial distribution of the LUE GPP values at the four percentiles over the globe in March–May, and those in June–August, September–November and December–February are provided in Figures S4–S6. It can be seen that there are obvious differences in the level of global GPP damage during different seasons. At the same time, three drought scenarios were set based on SPEI, moderate (−1.5 < SPEI ≤ 1), severe (−2 < SPEI ≤ −1.5) and extreme (SPEI ≤ −2). In this way, the temporal and spatial distribution characteristics of the probability of different degrees of damage to vegetation productivity under different drought scenarios were studied to determine areas around the globe susceptible to drought during different seasons.

As Figure 8d,h,i shows, the proportion of high-probability areas in the global terrestrial ecosystem increases with the increase in drought intensity when the GPP is lower than the 40th percentile in March–May. Under moderate drought conditions, the area with a probability of occurrence higher than 75% covers only 5.04% of the global terrestrial vegetation area, while under extreme drought conditions, the proportion increases to 27.78%. Further, when the GPP is lower than the 10th, 20th, and 30th percentiles, it still shows the same characteristics of change (the area with a probability of occurrence higher than 75% increases from 0.005% to 2.13%, from 0.12% to 9.24%, and from 1.65% to 18.71%). In other seasons, as the degree of drought increases, the conditional probabilities that GPP will be damaged at different levels increase as well (Figures S7–S9). Under the same drought severity with different levels of GPP damage, drought’s influence on GPP loss probabilities weaken gradually as the GPP damage level increases (Figure 8a–d). We explored the seasonal evolution in vegetation productivity’s loss probability further under the same drought scenario and the same level of damage of GPP. Using Figure 8h, Figures S7h, S8h and S9h as examples, it can be seen that under severe drought conditions, the area with a conditional probability higher than 75% was 17.48, 21.89, 17.37, and 14.25% in March–May, June–August, September–November, and December–February respectively, in which vegetation productivity was more sensitive to drought in June–August. The primary reason for this phenomenon is drought’s effect on GPP in mid- and high-latitude regions of the Northern hemisphere, as this period is the growing season of vegetation in those regions. The occurrence of drought during the growing season reduced the carbon sequestration capacity of vegetation greatly, and increased the GPP loss probability thereby.

Further, we analyzed the conditional probabilities of vegetation productivity losses under different drought scenarios in various climatic regions. Take the conditional probability when GPP is less than the 40th percentile from March to May as an example (as shown in Figure 9). Under the moderate drought scenario, the conditional probabilities of GPP loss are 0.54, 0.49, 0.45, and 0.44 in arid, semi-arid, semi-humid, and humid regions, respectively. (They are 0.74, 0.67, 0.56, and 0.56 under the severe drought scenario, and are 0.81, 0.72, 0.59, and 0.58 under the extreme drought scenario). Our results showed that arid and semiarid areas are more likely to suffer different levels of damage than are humid and semi-humid areas under different drought scenarios. As the degree of drought increases, the conditional probabilities of GPP loss also increase in different climate zones (similar results were also evident in other seasons (Figures S10–S12). By comparing the probability of vegetation productivity losses under different drought scenarios in different land cover types when the GPP ≤ 40th percentile (Figure 10, Figures S13–S15), it can be seen that as the degree of drought increases, the productivity loss probability of EBF, DNF, OS, SAV, and GRA show an increasing trend in different seasons. From March to August, for other land cover types, such as ENF, compared with extreme drought, severe drought has the greatest effect on the degree of damage to vegetation productivity, indicating that different land types respond differently to drought during different seasons.

4. Discussion

4.1. Estimating GPP and Drought’s Effects on GPP

In this study, special attention was given to determining water stress (including atmosphere vapor pressure deficit, soil moisture content and the humidity deficit) in the LUE model to estimate the global GPP. We found that different LUE models have a good fit-in estimating the LUE GPP (R² exceed 0.7). Cai et al. [42] showed that the range of global GPP estimated by different light energy utilization efficiency models was 95–140 Pg C yr⁻¹, which is consistent with this study’s estimation results. Generally, because of the lack of water during drought periods, vegetation’s stomatal conductance tends to decrease to minimize water loss and prevent hydraulic conductivity loss [43,44], which reduces GPP. Many studies have shown that drought reduces terrestrial ecosystem carbon sinks significantly, and may even transform them into carbon sources [13,45]. However, our results showed that approximately 30% of the global GPP values increased after being disturbed by drought (Figure 2b), a phenomenon attributable to climatic factors’ influence. Our results showed that GPP is correlated negatively with temperature and solar radiation in arid and semiarid areas, while the converse is true in humid and semi-humid areas. Further, GPP is correlated positively with soil moisture globally. When a drought episode occurs, because of lower radiation and temperature, and the replenishment of soil moisture in the high latitudes of the Northern Hemisphere, Eastern Asia, the Amazon River basin and Central Africa, it causes an increase in GPP in these regions relative to those in non-dry periods (as shown in Figure 11). This result can be confirmed by Zhao and Runing [46]. In contrast, in most parts of the world, GPP decrease during drought periods because of the increase in radiation and temperature and the decrease in soil moisture.

4.2. Terrestrial Ecosystems’ Drought Resistance

Terrestrial ecosystems’ drought resistance can be characterized as: a water deficit must continue for a period of time before negative anomalies occur in ecosystem variables [47]. In this study, response time of productivity to water availability was defined according to timescales over which the maximum GPP-SPEI correlations were recorded. To some extent, the response time of vegetation productivity to drought can characterize terrestrial ecosystems’ drought resistance. Our research results showed that 56.8% of the global terrestrial ecosystems’ response time to water availability is based largely on short- and medium-term time scales. Because different hydrological backgrounds and vegetation composition affect terrestrial ecosystems, their drought resistance demonstrates strong spatial heterogeneity. In Canada, Northwestern Colombia, Eastern Russia, and Northeastern and Southeastern regions China, the response time of vegetation productivity to SPEI is relatively long, usually greater than 12 months. The longer the response time to water availability changes in these regions, the stronger the vegetation productivity’s capacity to withstand long-term water shortages. Seasonal changes also affect the resistance of the ecosystem to drought in the same region, particularly in mid-to-high latitude regions (Figure S3). Taking the mid-to-high latitude regions of the Northern hemisphere as an example, it can be seen that the proportion of GPPs in the high latitude with a longer response time to SPEI decreases gradually with the passage of seasons. This result indicates that vegetation has a strong resistance to drought in the early stage of growth, and as subsequent vegetation growth requires considerable water, the occurrence of drought at this stage is more likely to cause vegetation productivity to decline.

Climatic conditions and different types of land cover often affect terrestrial ecosystems’ drought resistance [48,49]. Our results showed that as the climatic conditions become gradually humid, the response time of vegetation productivity to changes in water availability increases. This indicates that vegetation productivity in arid and semiarid areas is more vulnerable to drought than in humid and semi-humid areas. Vicente-Serrano et al. [49] also found that arid biomes respond to drought at short time-scales. This may be related to the adaptation strategies of vegetation to water availability in different climate zones. When a water deficit occurs, arid ecosystems adapt quickly to water shortages by reducing water loss and respiratory costs and increasing growth rates. Fang et al. [20] demonstrated that, as there is a long-term water surplus in humid areas, a short-term drought can not cause changes in vegetation productivity easily. Further, our results revealed different land biomes’ role in drought resistance. By comparing GPP’s response time to SPEI under different land cover types, we found that DNF, OS, EBF, and ENF have a greater ability to adapt to drought, while WSA, MF, SAV, and GRA have a weak ability to resist and mitigate drought. Forests generally have deep root systems, and thus, under severe drought conditions, they can use the water stored in the deep soil to achieve strong resistance to drought. Whether drought disturbs a forest ecosystem is related closely to its intensity and duration. The water use strategy of shrubs is generally to increase the use of surface soil water during non-arid periods and absorb deep soil water during dry periods to maintain their own production needs [50]. However, herbaceous plants’ xylem system has low water and carbon storage capacity compared with that of woodland and shrubs, so it is less resistant to drought. Ivits et al. [51] found that steppic ecosystems also showed weak resistance against drought. Our results are more consistent with the research above.

4.3. The Significance for Ecosystem Management of Estimating the GPP Loss Probability

We developed an optimal bivariate probabilistic model to derive the vegetation productivity loss probabilities under different drought scenarios using copula method. The vegetation productivity loss probability drought causes can also reflect the ability of terrestrial ecosystems to resist interference. Under the same drought conditions, the higher the probability of vegetation loss, the weaker the resistance to drought. Our results showed that arid and semiarid areas have higher conditional probabilities of vegetation productivity losses under different drought scenarios. The productivity loss probability of EBF, DNF, OS, SAV, and GRA showed an increasing trend during different seasons. Quantifying the probability of varying degrees of damage to vegetation productivity under predictable drought scenarios has great significance in mitigating and adapting to global changes. However, there is none of the comparisons either between vegetation are types or the prediction is self are not statistically significant. Huang et al. [52] found that the persistence of a water deficit (11 months) with an intensity of −1.64 (SPEI) led to negligible growth of conifer species. Berdugo et al. [53] showed that aridification can lead to systemic and abrupt changes in multiple ecosystem attributes. Therefore, in future studies, to determine how much drought can cause significant changes in ecosystem productivity and to estimate the corresponding GPP loss probability can better reflect the differences between arid to humid areas.

5. Conclusions

In this study, we explored satellite-observed global terrestrial vegetation production in response to water avaliability by determining gobal vegetation productivity’s seasonal dynamics in response to drought in various climate zones and land biomes and quantifying its vulnerability under predictable drought scenarios. Our primary conclusions are as follows:

• Different LUE models have a good fit effect in estimating GPP. The fitting R² of VPD_GLO-SM, VPD_MOD-SM, VPD_GLO-ETR and VPD_MOD-ETR were 0.7739, 0.7399, 0.7427, 0.7459 and 0.7628, respectively. From 1982 to 2015, the global mean annual GPP of terrestrial vegetation continued to increase at an average rate of 0.134 Pg C a ⁻¹ (p < 0.001), but its growth rate declined after the mid-1990s. GPP is expected to decrease in 71.91% of the global land vegetation area because of increases in radiation and temperature and decreases in soil moisture during drought periods.
• Vegetation productivity and water availability are largely correlated positively globally. Further, seasonal changes also affect vegetation productivity’s dependence upon water availability. The correlation coefficient between GPP and SPEI declined from 0.76 to 0.47 as the climatic conditions became gradually humid, indicating that the vegetation productivity in arid and semiarid areas depends more heavily on water availability than that in humid and semi-humid areas. Various land cover types have different adaptation strategies to the increase and loss of water resources, and the productivity of GRA, SAV, and DBF has a higher correlation with water availability.
• 56.8% of the global terrestrial ecosystems’ response time to water resources is based primarily on short and medium-term time scales (3–6 months). The GPP’s mean response time to SPEI increased from 3.9 to 8.9 months as the climatic conditions became gradually humid, which indicates that the capacity of productivity of vegetation in arid and semiarid areas to withstand long-term water shortages is weaker than that in humid and semi-humid areas. The land cover types that are more relevant to water availability are often accompanied by weak drought resistance, while DNF, OS, EBF and WET have a stronger ability to resist long-term water deficits.
• Under the scenario of the same level of GPP damage with different drought degrees, as droughts increase in severity, GPP loss probabilities increase as well. Further, under the same drought severity with different levels of GPP damage, drought’s effect on GPP loss probabilities weakens gradually as the GPP damage level increases. Similar patterns were observed in different seasons. Our results showed that arid and semiarid areas have higher conditional probabilities of vegetation productivity losses under different drought scenarios. The productivity loss probability of EBF, DNF, OS, SAV, and GRA show an increasing trend in different seasons, and different land types have different responses to drought in different seasons.