Compound drought constrains gross primary productivity in Chinese grasslands

Water constraints disturb and damage the growth and development of grassland vegetation mainly through both atmospheric and soil pathways. In the background of rapid climate change in the future, the impacts of water constraints on grasslands are expected to further deepen. However, current studies lack reports exploring the frequency, intensity, and area of land-atmospheric compound drought on carbon indicators in grassland ecosystems. In this study, we analyze the response of China grasslands to dual terrestrial-atmospheric water constraint events using ISIMIP gross primary productivity (GPP) data to reveal the carbon cycle-climate feedback relationships over the Chinese grassland. We found that the occurrence probability of compound drought events (i.e. land-atmospheric water constraint) was 3–4 times higher than that of random drought events, and the frequency, intensity, and affected area of compound droughts were significantly higher than that of single droughts. Compound droughts caused a decline of up to 20.27% in GPP of grassland ecosystems in China, while the decline of single atmospheric drought or soil drought was only 12.34% and 14.32%. Which is because vapor pressure deficit and soil moisture are a set of strongly coupled bivariate variables, and the continued strengthening of the land-atmospheric feedback causes a higher probability of occurrence of compound drought events and an increased impact on ecosystem GPP.


Introduction
Soil drought and atmospheric drought are considered to be two hydraulic processes that affect ecosystem productivity during drought (Zhou et al 2019a). Plants initially draw water molecules from the soil to provide water for photosynthesis, while simultaneously dispersing them in the atmosphere through transpiration of water molecules from plant leaves (Shen et al 2015). Both of these stresses significantly reduce terrestrial carbon uptake and crop yields, and cause widespread vegetation mortality (Zhang et al 2016, Rigden et al 2020. Soil moisture (SM) is the most important source of water directly available to plants and is an important monitoring indicator reflecting drought (Li et al 2018). Inversion and simulation of SM information through satellite remote sensing, land surface models and SM balance models have been developed considerably (Li et al 2018, Sehgal et al 2021, and can provide direct data sources for large-scale monitoring of the surplus and deficit of water required by vegetation. The effects of reduced SM on terrestrial vegetation are now well established, with SM shortage suppressing water supply to vegetation (Lian et al 2020, Liu et al 2020, and being a direct factor leading to reduced photosynthetic capacity or even death of vegetation (Liu et al 2019). Meanwhile, when terrestrial ecosystems face long-term or short-term droughts, plants can excessively deplete water in the surface and shallow soils (Ciais et al 2005), exacerbating the intensity and probability of occurrence of soil drought.
Meanwhile, when atmospheric water is deficient, plants tend to inhibit photosynthesis by reducing leaf stomatal conductance and may exacerbate soil drought by enhancing evapotranspiration, exposing vegetation growth to severe drought stress (Zhang et al 2014, Sinclair et al 2017. However, there is little monitoring and assessment of atmospheric drought, because conventional monitoring tools and approaches have difficulty in accurately capturing the dynamics of atmospheric water (all water evaporated from the ocean and land) (Lobell et al 2014). The vapor pressure deficit (VPD) provides an approximate description of the atmospheric moisture gain or loss as the difference between the pressure exerted by water vapor when the air is fully saturated and the actual pressure exerted (López et al 2021). When the difference increases, the stomatal conductance of the leaves is reduced to prevent water dissipation in the plant, but this is at the expense of reducing the photosynthetic rate (Grossiord et al 2020). A recent report stated that global vegetation was subject to atmospheric drought that has been intensifying since 1990, resulting in an ecosystem whose carbon uptake capacity has been greatly reduced (Yuan et al 2019).
Here, from the two water potential processes of water uptake and water consumption by plants, a drying event that occurs simultaneously in the soil and atmosphere at the same place at the same time in an ecosystem is defined as a compound drought. Compound drought is driven by a series of complementary land-atmosphere physical processes (Zhou et al 2019a). On the one hand, low SM reduces terrestrial evapotranspiration, leading to higher temperatures and VPD, which are associated with near-surface desiccation due to evaporative cooling. On the other hand, elevated VPD promotes enhanced terrestrial evapotranspiration, which in turn forces SM dissipation (Massmann et al 2019). With global warming, extreme land-atmosphere compound droughts occur frequently and have a great impact on ecosystem structure and function.
Currently, studies on the relative impacts of landatmospheric water constraints on ecosystem gross primary productivity (GPP) are scarce and have focused singularly on analyzing the effects of SM (or atmospheric moisture) surpluses and deficits on vegetation. In particular, it is rare to assess the probability of simultaneous occurrence of two water constraints (atmospheric drought and soil drought), which are most critical for plant growth in Chinese grasslands.
Therefore, it is important to enhance understanding of water stress effects of land-atmospheric compound drought in Chinese grassland ecosystems, to point out the probability that soil water deficit and atmospheric water deficit events may occur together, and to explore the relative intensity and extent of land-atmospheric water constraints on GPP in Chinese grassland ecosystems, which is an important theoretical reference for drought disaster emergency prevention and policy formulation in China. Based on this, the research objectives of this paper as follow: (a) Using correlation and Granger causality tests to explore the dependence structure between SM and VPD. (b) A probabilistic evaluation framework was developed to assess the probability of simultaneous occurrence of land-atmospheric water constraint events by introducing couplas functions. (c) To analyze the impact of the area, frequency, and intensity of the occurrence of compound droughts on the Chinese grasslands GPP compared with that of single droughts.

Model GPP products, simulation settings
We used four global vegetation dynamic models (CARAIB, LPJ GUESS, LPJmL, ORCHIDEE DGVM) to output the results: GPP (monthly value, 0.5 • resolution). The model GPP products were all based on the uniform atmospheric circulation model IPSL-CM5A-LR as meteorological forcing data, and considered two scenarios: historical period (including the impact of human emissions, 1901-2005), medium-high greenhouse gas (GHG) emission scenario (RCP6.0, 2006(RCP6.0, -2099. The GPP simulation results of all models are performed strictly under the ISIMIP2b standard protocol (www.isimip.org/), so the difference between model output results is only related to the complexity and use of the model. Considering that a single model provides valuable and usable insights only at the regional scale, to effectively eliminate intra-model variability and reduce uncertainty in GPP products, we do pooled averaging of multi-model GPPs for the next step of analysis.

SM products
SM data provided using the ISIMIP2b simulation protocol was used to characterize SM wetness and dryness. Our selected the set of average SM products output by four models. However, that different models differ in simulating land-atmospheric exchange fluxes and carbon-water cycle reserves in natural and agro-ecosystems, so each model provides a different thickness of SM. To match the soil thickness and reduce the data error, the SM we calculated was limited to the sum of the water fractions of all soil layers within 3 m underground.

VPD calculation
We also obtained standardly corrected model input parameters (temperature and relative humidity) from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) for calculating VPD. VPD refers to the difference between saturated vapor pressure (SVP) and actual vapor pressure at a specific temperature and is a direct measure of the intensity of atmospheric drought (Buck 1981) SVP is a nonlinear function of air temperature, which can be obtained directly from the calculation of temperature, and the empirical formula is as follows: Ta+243.5 , (2) wherein, T a is the surface temperature of the air ( • C); Z is the elevation (m); P msl is the atmospheric pressure at mean sea level (≈1013.25 hPa); P mst is the atmospheric pressure (hPa)

Remote sensing GPP products
Modeled GPP is not independent of SM and VPD because they are already 'simulated' with inherent functions of SM and VPD. The GPP products of multiple models provided by ISIMIP are not independent of the SM and VPD data. Therefore, we have tried to validate the conclusions obtained from the model GPP using the GPP dataset from satellite observations in order to avoid the over-simulation of the model GPP. Here, we used two of the more novel and mainstream GPP products.  Sun-induced chlorophyll fluorescence (SIF) is a new remote sensing technique that has been rapidly developed in recent years, and the close connection between SIF and photosynthetic processes makes it an effective probe for indicating vegetation photosynthetic changes and a powerful tool for monitoring GPP (Wang et al 2020b). A novel vegetation index (near infrared reflectance of terrestrial vegetation, NIRv) obtained based on remote sensing data, i.e. the product of normalized difference vegetation index (NDVI), a normalized vegetation index, and NIR band reflectance, is highly correlated with remotely sensed SIF products; the results based on mechanistic derivation, model simulation, and analysis of remote sensing data all show that NIRv can be used as an alternative to SIF for estimating global GPP. Therefore, in the validation, previous studies used NIRv-based GPP product (0.05 • , monthly) provided by Wang et al (2020b) for validating the reliability of the model GPP. GPP NIRv is available from the URL: http://data.tpdc.accn/zh-hans/data/d6dff40f-5dbd-4f2d-ac96-55827ab93cc5/?q=gpp to download.

Global ecosystem GPP long time series data (1982-2018) based on remotely sensed EVI and light use efficiency model
The Global Land Surface Satellite (GLASS), the first set of advanced global land remote sensing products with independent intellectual property rights in China, was widely used in the study of global change. Among them, the vegetation GPP is the most important product dataset of GLASS, which spans nearly 40 years  and is the primary choice for studying the long-term evolution characteristics of GPP due to its relatively high spatial and temporal resolution (0.05 • , monthly). The earlier model was driven by only four variables: NDVI, photosynthetically active radiation, air temperature, and the Bowen ratio of sensible to latent heat fluxes (used to calculate water stress) (Yuan et al 2007). It was noted that GLASS GPP was shown to better characterize photosynthetic productivity of major global ecosystems compared to other satellite-based GPP models (Yuan et al 2014). The link to access this data is: www.glass. umd.edu/Download.html.

Screening of GPP, SM and VPD in the warm season
Previous defined the warm season in grid cells as the average of the three hottest months (one value per year) or the three months with the highest mean temperature (three value per year), given that the warm season coincides with the main growing season of plants (Zscheischler andSeneviratne 2017, Zhou et al 2019b). We will analyze the average (one value per year) GPP and VPD dependence for mean value of the three months with hotter temperatures, as carbon loss due to moisture shortage tends to be most dramatic in the hottest three months. Therefore, the locations of the three months with the highest temperatures on the grid cell were extracted from the monthly mean temperatures for the period 1901-2005, and the GPP, SM, and VPD for the three months with the highest temperatures on the pixel scale for the historical and future periods were selected in the time series for further analysis.

Interannual correlation measures
For bivariate extremes, correlation is a good indicator and is commonly used to assess the possible impact of structural dependence on binary extremes (Zscheischler and Seneviratne 2017). We used an model dataset for calculating interannual correlations between VPD and SM. In addition, to remove the influence of climate change signals on long-term trends, we performed linear detrending of the bivariate prior to calculating correlations. In this paper, the Spearman rank correlation coefficient was chosen to analyze the correlation between VPD and SM, with a stronger negative correlation of the bivariate indicating more significant negative feedback.

Granger causality test
A one-sided correlation analysis is not comprehensive enough to assess the bivariate coupling, so we then introduce the Granger causality test to assess the dependency structure between the VPD and SM bivariates. The Granger causality test was originally used by economists to analyze the causal relationship between two variables (Granger 1969). Because it takes into account the variation in the variables themselves when analyzing the relationship between variables, it is good at removing pseudo-correlation between variables, and the method has been widely used in ecological causality studies in recent years (Xie et al 2019). The Granger causality for a long time series of bivariate variables X, Y can be described as follows: when the data contains information on the history of variables X, Y, the prediction result for variable Y is better than the prediction result for Y considering only the historical information from Y. In other words, when information on variable X can contribute to explaining future changes in variable Y, then variable X is considered to be the Granger cause of variable Y. It follows that the Granger causality test is essentially a prediction of the causal link between two sets of smooth time series data. The vector autocorrelation regression model is described as follows where P is the maximum lag order; N is the sample size; α and β are the regression coefficients; ε 1 and ε 2 are the error terms, and both are assumed to be uncorrelated.
The constrained F-test is used to test whether the Granger test is passed RSS m and RSS r are the sum of squared residuals from the regression model and F-test respectively. The original hypothesis that X is the Granger cause of Y is not rejected if F ⩽ F α . where F α is the confidence level (set at 0.05).
We regress the Granger causality between VPD and SM on an pixel-by-pixel basis. If using VPD as a predictor improves the predictive accuracy of SM, and GPP as a predictor but does not improve the predictive accuracy of VPD, then VPD is considered to be a one-way Granger cause of SM. Conversely, the SM is considered to be a one-way Granger cause of the VPD. If VPD is the Granger cause of SM, and SM is the Granger cause of VPD, then there is a two-way causal relationship between VPD and SM.

Bivariate linkage to calculate the probability under compound drought
We constructed a probabilistic framework for assessing the probability of simultaneous occurrence of atmospheric and soil drought using the Couplas function based on the method and steps proposed by Zscheischler and Seneviratne (2017) as follows.

Marginal distribution fit
Bivariate frequency analysis requires determining the distribution of random variables U and V, so determining the marginal distribution of bivariates is a prerequisite for constructing a joint probability distribution. Kernel density estimation is the most widely used test method in the field of nonparametric estimation (Jones and Hardle 1992), where a kernel distribution produces a nonparametric probability density estimate that adapts itself to the data, rather than choosing a probability density estimate with specific parameters. We used kernel density estimation to find marginal distribution fits for VPD and GPP.
The calculation formula of the kernel function density estimation method is described as follows: f n (x) represents the kernel density value; n represents the number of samples in the bandwidth range; h is the window or bandwidth, represents a reasonable smoothing parameter, and h > 0; K is the kernel smoothing function; (x-x i ) is the distance from point to x i , x 1 , x 2 , …, x n is a random sample from an unknown distribution. The kernel density estimate of the cumulative distribution function is described as follows:

Fitting and preferring of joint probability distributions
Copulas is a multivariate distribution function with a domain of [0,1] that is used to describe the correlation between multiple variables. The binary variable copulas function is often used to describe the dependency structure between two sets of random variables and to count the joint probability of an event (e.g. a compound drought) occurring. Sklar confirms that the copula function is unique in the sense that if F(·,·) is a joint distribution function and X(·) and Y(·) are marginal distribution functions of independent variables, then there must exist a copula function C(·,·) that satisfies When X(·) and Y(·) When continuous, there must be a uniquely deterministic C. Conversely, X(·) and Y(·) are only one-dimensional distribution functions.
According to Sklar's theorem, the joint probability distribution function F X,Y (x, y) for variables X and Y can be expressed as where Fx(x) = P(X ⩽ x) and Fy(y) = P(Y ⩽ y) are the cumulative distribution functions of variables X (e.g. VPD) and Y (e.g. −GPP). C is the joint distribution function of U = Fx(x) and V = Fy(y) after marginal fitting, and the new sequences U and V after the marginal distribution fit transformation have the characteristics of uniform distribution. For example, the probability that both variables of an event exceed a given threshold can be expressed as: The copula family junction function is used to reconstruct the dependency structures of VPD and −SM. To more accurately describe the dependency structure of binary variables, the fitted optimal one is selected from five joint distribution functions of Clayton, Frank, Gumbel, t, and Gassian Copula for further analysis. The method of the goodness test is based on the minimized squared euclidean distance, which is described as follows: where CUV is the empirical value of each binary copula function fitted, and C is the theoretical value.

Land-atmospheric compound drought transcendence probability calculation
We define a compound extreme event as one in which the VPD is above the 90th percentile of the time series while the SM is below the 10th percentile of the time series (equivalent to the −SM being above its 90th percentile). The joint probability of this compound extreme event is as follows p = P(U > 0.9 ∩ V > 0.9) = 1 − F (0.9) − F (0.9) + C (0.9, 0.9) = 1 − P(U > 0.9) − P(V > 0.9) Based on equations derived (Serinaldi 2014, Guo et al 2020) the probability of simultaneous occurrence of atmospheric drought and soil drought (complex drought) is For a bivariate with a 10% extreme threshold, the probability of a random sample falling in a 10 × 10 grid is 1%. The initial probability of the bivariate is considered to be exceeded if the probability of both atmospheric and soil drought occurring at the same time exceeds 1%, which is defined as the exceedance probability in this paper.

Dependency structure of SM and VPD bivariate
Land-atmospheric feedbacks are potentially linked to the occurrence of both soil drought and atmospheric drought, and it is critical to explore the role of landatmospheric feedbacks in the occurrence of compound drought events. After analyzing the bivariate dependence between SM and VPD in the warm season, we found that SM and VPD in Chinese grasslands showed negative correlations at meta-scale for both historical and future periods, with correlation coefficients of −0.364 and −0.547, respectively (figure 1). In terms of region, the bivariate correlation of the Loess Plateau is the most significant with negative correlations of −0.411 and −0.612 for historical and future periods, respectively, followed by Inner Mongolia and Xinjiang, and negative correlations of −0.366 and −0.500 for the Qinghai-Tibet Plateau, respectively. This confirms that land-atmospheric is indeed a strongly coupled set of variables, and reflects that the correlation between SM and VPD is a good indicator to explore the possibility of extreme soil drought and atmospheric drought simultaneously.
We then used Granger causality analysis to detect potential feedbacks between SM and VPD during the warm season, significant bivariate causality was found in more than 65% of areas in both historical and future periods in Chinese grasslands (figure 2), indicating detectable mutual feedbacks between the SM and VPD bivariates. VPD was considered to be a unidirectional cause of SM in more than 25% of the areas, and these pixels were mainly scattered across Chinese grasslands. The unidirectional causes of SM  The red area is considered as a two-way causality between GPP and SM; green indicates that VPD is a one-way cause of SM; and blue indicates that SM is a one-way cause of VPD. as a cause of VPD occurred locally in Inner Mongolia and Xinjiang, accounting for less than 5% of Chinese grasslands. In conclusion, above results shown that there is a strong feedback relationship between landatmospheric in Chinese grasslands, and that atmospheric drought triggers secondary soil drought in a wider area.

Frequency, intensity and impact of compound drought events on ecosystem GPP
The results in section 3.1 indicated that the probability of simultaneous occurrence of atmospheric drought and soil drought is high, but the multiplier of the probability of occurrence of compound drought is unknown. Considering that the couplas function can better associate the tail correlation of quantitative extreme variables , Wang et al 2021, the likelihood multiplication factors (LMFs) of simultaneous occurrence of atmospheric droughts and soil droughts based on the couplas probability framework were shown in figure 3. The results shown that the probability of occurrence of compound drought in Chinese grasslands is much higher than the probability of occurrence of two events at random (0.1 × 0.1 = 0.01, i.e. 1%, for a random event of two variables). The LMF for Chinese grasslands in the historical period averaged 3.085% and 4.234%, respectively (figures 3(a) and (b)), implying that the probability of a compound drought occurring simultaneously is more than 3 ∼ 4 times higher than the probability of a random event. The mean value of the difference between the bivariate correlations of VPD and SM for the future and historical periods for the warm season is −0.182% (figure 3(c)), and the difference in the LMF is 1.149% ( figure 3(d)), which indicates that the probability of extreme compound drought events is significantly higher (over 1.149%) in the future than in the historical period. In summary, above result reflects the much higher frequency of simultaneous land-atmospheric moisture constraints compared to the scenario assumed by the independent combination of univariate statistics, which can have direct and far-reaching impacts on ecosystems.
When compared the relative anomalies of composite drought, atmospheric drought, and soil drought leading to the decline of ecosystem GPP in Chinese grassland. We found that the magnitude of GPP decline caused by compound drought is higher than that of single drought, and compound drought can lead to a decline of 20.27% and 34.04% in ecosystem GPP in historical and future periods, respectively (figure 4). Atmospheric droughts caused lower decreases in GPP, with 12.34% and 23.08% decreases, respectively. The decline in ecosystem GPP under soil drought constraint was slightly higher than that of atmospheric drought, with 14.32% and 25.18% declines in historical and future periods, respectively. The analysis revealed that the threat of compound drought on carbon uptake and carbon accumulation in Chinese grassland ecosystems is clearly higher and more influential compared to single drought.

Similarties and differences of combined drought on model GPP and satellite GPP
We mentioning that the remote sensing-based GPP products are limited by the age of satellite observations, and we can only obtain datasets since 1982. Moreover, the spatial resolution of the remotely sensed GPP dataset is much higher than that of the model GPP dataset. Therefore, before conducting the validation formal analysis, we preprocessed the data. First, we recalculated the bivariate correlations and probability multiplication factors for VPD and SM from 1982 to 2018 to match the deficiency of the remotely sensed dataset that only provides nearly four decades of observations. Considered that the ISIMIP dataset is divided into historical and future periods with the 2005 time node, we integrated the VPD and SM for the years of historical period  and future period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) in the same time series before the calculation. Second, because the remote sensing data resolution is 0.05 • , we resampled the GPP NIRv and GPP GLASS datasets with 0.05 • resolution to 0.5 • by the nearest neighbor method before performing the validation.
Our results found a strong negative correlation between the bivariate VPD and SM from 1982 to 2018, and a very high probability of extreme drying events occurring simultaneously in VPD and SM (figure 5). This suggests that both atmospheric drying and soil drying may have strongly weakened the ecosystem GPP during this time period.
To verify the possibility of model GPP, figure 6 showed the responses of GPP ISIMIP , GPP NIRv , and GPP GLASS to compound drought, atmospheric drought, and soil drought, respectively, during 1982-2018. The results showed that compound drought presents a strong constraint on the GPP of Chinese grassland ecosystems regardless of model GPP or remotely sensed GPP, and the effect of this moisture constraint, compound drought events, far exceeds that of a single atmospheric drought or soil drought. This is extremely similar to the thematic conclusions presented in this study.
However, peculiarly, we also found that the two remote sensing-based GPP products were much less sensitive to drought event responses than the model GPP. The remote sensing GPP did not seem to present   well enough to produce very strong feedbacks to drought event responses in the Chinese grassland GPP.

Discussion
Climate anomalies have fundamentally changed the inherent variability of weather phenomena such as precipitation and temperature, leading to the frequent occurrence of extreme weather such as droughts and heat waves (Ciais et al 2005). In particular, land-atmospheric compound droughts have become the most concerned natural disasters because of their highest frequency and impact Running 2010, Zhou et al 2019b). Previous studies have often considered the effects of high VPD and low SM events on ecosystem GPP to be independent, but in fact high VPD and low SM events often occur simultaneously. We found that the probability of simultaneous atmospheric and soil droughts in Chinese grasslands was 3 ∼ 4 times higher than that of stochastic droughts, implying that land-atmospheric compound moisture constraint events occur more frequently and are more extreme in intensity of exposure ( figure 3(a)). Moreover, high VPD and low SM events are expected to occur more frequently and more extremely in the future ( figure 3(b)). The reason for the apparent increase in the intensity and frequency of compound droughts may be due to the strong feedback between land-atmospheric that tends to be tighter, and this strong feedback is bound to be stronger in the future. Granger's causality test reveals a bidirectional causal relationship between VPD and SM in most areas of Chinese grasslands, which may be due to elevated VPD actively driving terrestrial evapotranspiration to accelerate SM depletion, while reduced SM decreases evapotranspiration to indirectly raise temperature and VPD (Berg et al 2016, Green et al 2017, Zhou et al 2019a, the strong feedback between land-atmospheric drives the dramatic outbreak of compound drought events. Zhou et al (2019a) experiments using GRACE-CMIP5 showed that simultaneous soil drought and atmospheric dryness are greatly exacerbated by land-atmosphere feedback. The feedback of soil drought to the atmosphere has largely contributed to the extremeization of atmospheric dryness. In addition, SM-precipitation feedback amplifies precipitation and SM losses in most areas. Therefore, we state that it is the strong coupling of VPD and SM that results in a high probability of soil drought when atmospheric drought occurs, and a correspondingly high probability of secondary atmospheric drought triggered by soil drought.
When compared ecosystem responses to extreme events at the meta-scale, positive or negative GPP anomalies during compound or independent droughts were analyzed. We found that compound extreme VPD and SM reduced terrestrial carbon uptake more strongly than single drought stresses considered alone in terms of GPP as a carbon indicator, i.e. compound drought events had an increasingly negative impact on ecosystem GPP (figure 4). This may be due to the fact that moisture status is the largest limiting factor affecting global vegetation growth and carbon cycling (Xiao et al 2013), and the superimposed drought from the atmosphere and soils multiplies the moisture stress, more strongly and directly reducing the availability of water required by vegetation (Zhang et al 2016). In the future, it is foreseeable that the frequency and intensity of the compound extreme events will further increase as temperatures continue to increase (Wang and Chen 2014). Crucially, increasing extreme compound moisture constraints may offset the effects of carbon and nitrogen fertilization (Reich et al 2014, Wang et al 2020a, and directly threaten the ability of grassland ecosystems to act as important carbon reservoirs and sinks. We, therefore, recommend that governments and international organizations adopt and develop stronger policies to mitigate the potentially widespread and profound impacts of compound drought on terrestrial ecosystems to the extent possible. Upon analyzing the response of Chinese grassland ecosystems to desiccation events in the historical and future periods, we found that very high VPD and very low SM become more extreme in the future due to the concomitant increase of greenhouse gases, i.e. compound extreme events become intensified in the future. The LMF for simultaneous landatmospheric compound drought events in the future period was 4.149%, higher than that of 3.02% in the historical period (figures 3(a) and (b)). The percentage of area with LMF > 5% was only 4.18% in the historical period, but increased to 25.87% in the future period. When compared the probability differences between the two periods, 62.54% of the areas had LMF increases greater than 1%, and 12.45% of the areas had LMF increases greater than 3%. In particular, the LMF variation of grasslands on the Loess Plateau and Xinjiang was greater ( figure 3(d)). Obviously, the expansion of the area and the increase of the frequency of the future compound drought range pose significant challenges to the ability of each region of Chinese grasslands to act as a carbon sink in the future. From figure 4, we found that the intensity of the decline in GPP of Chinese grasslands due to compound drought in the future period can be higher than that of the historical period by about 13.77%. In conclusion, the area, frequency, and intensity of compound droughts in Chinese grasslands in the future period are higher than those in the historical period, and it is expected that the increase in the intensity and frequency of compound drought events may correspond to an increased risk of plant mortality, which is consistent with the results of previous studies (Zhou et al 2019b).
After analyzing the impact of compound drought on grassland productivity in various regions, we found that the grassland productivity of the Qinghai-Tibet Plateau in the historical period showed a different change after undergoing dual water stress of land and atmosphere. Specifically, grasslands in Inner Mongolia, Xinjiang, and the Loess Plateau showed severe GPP losses after compound droughts, but grassland GPP on the Qinghai-Tibet Plateau showed a certain degree of increase ( figure 4(a)). We believed that the Tibetan Plateau has a unique alpine ecosystem, which may be the reason why the conclusions differ from other regions. Over the past 100 years, the rate of warming on the Tibetan Plateau has been much higher than in other regions (Zhang et al 2021), which has significantly contributed to the upward GPP on the Tibetan Plateau grasslands in both ways. First, warming will significantly lift the temperature limits of alpine ecosystems, which will greatly release the potential for plant growth, which in turn will promote GPP increase. Second, the warming forces the otherwise solid snow-ice to melt, which allows the soil of the Tibetan Plateau to seep in large amounts of water needed (Deng et al 2019). Although the frequency and intensity of compound water constraints have increased, it is still difficult to offset the GPP increase caused by the warming of plateau ecosystems. However, we must point out that the GPP increase brought about by this long-term warming is not sustainable. We also found that in the coming period, the Tibetan Plateau, like other grasslands in China, will have compound droughts that limit the decline in GPP ( figure 4(b)). This is mainly due to the fact that the effect of warming to stimulate plant growth is short-lived, and when the temperature reaches a certain level in the plateau ecosystem, the grassland productivity associated with the warming will be converted from promotion to inhibition. Moreover, the total amount of ice-snow on the Qinghai-Tibet Plateau is fixed, and when the icesnow disappear, the water source that plants can draw from the soil will also decline rapidly, which will undoubtedly inhibit the accumulation of productivity in the future.
Furthermore, discerning the effects of compound droughts on ecosystems remains challenging, and even a probabilistic perspective can only provide possibilities but not definitive conclusions. This is because the hydraulic response of plants to soil stress and atmospheric demand is highly simplified in models, which is a major source of uncertainty in terrestrial carbon response to climate change (Zaehle et al 2005). Further clarification of this long-standing and complex issue is needed to open new avenues for improved models and better management of compound drought risk.

Conclusions
Used four terrestrial ecosystem models, we found that (a) high VPD and low SM events are a set of strongly coupled bivariate variables, and enhanced land-atmospheric feedbacks may be responsible for more frequent than expected composite events of high VPD and low SM. (b) The impact of landatmospheric composite drought events on Chinese grassland GPP is much higher than the strong negative impact of a single drought on Chinese grassland productivity. Our study highlights the very broad threatening nature of compound extreme drought events on the ability of Chinese grasslands to act as carbon sinks.