Enhanced dominance of soil moisture stress on vegetation growth in Eurasian drylands

ABSTRACT Despite the mounting attention being paid to vegetation growth and their driving forces for water-limited ecosystems, the relative contributions of atmospheric and soil moisture dryness stress on vegetation growth are an ongoing debate. Here we comprehensively compare the impacts of high vapor pressure deficit (VPD) and low soil water content (SWC) on vegetation growth in Eurasian drylands during 1982–2014. The analysis indicates a gradual decoupling between atmospheric dryness and soil dryness over this period, as the former has expanded faster than the latter. Moreover, the VPD–SWC relation and VPD–greenness relation are both non-linear, while the SWC–greenness relation is near-linear. The loosened coupling between VPD and SWC, the non-linear correlations among VPD–SWC-greenness and the expanded area extent in which SWC acts as the dominant stress factor all provide compelling evidence that SWC is a more influential stressor than VPD on vegetation growth in Eurasian drylands. In addition, a set of 11 Earth system models projected a continuously growing constraint of SWC stress on vegetation growth towards 2100. Our results are vital to dryland ecosystems management and drought mitigation in Eurasia.


Vegetation growth indices
Three satellite-observed vegetation indices as a proxy for vegetation condition were used in this study, including the third-generation biweekly Advanced Very High-Resolution Radiometer (AVHRR) normalized difference vegetation index (NDVI) (GIMMIS-NDVI3g) [1,2], the daily Ku-band vegetation optical depth (VOD) [3], and gross primary production (GPP) dataset based on NIRv (GPP NIRv ) [4]. We used three independent remote sensing datasets because each has its own limitations. For instance, NDVI is affected by soil background; VOD is influenced by soil moisture; and GPP depends on meteorological factors, furthermore, the fitness between NIRv and GPP has large uncertainty. All vegetation data were resampled to a 1° spatial resolution (Supplementary Table 2). The bi-weekly GIMMIS-NDVI3g data was composited into monthly values by maximum value compositing method. We only used vegetation index values averaged over the growing season period (April to October) to avoid spurious vegetation changes caused by snow presence in dormant seasons [2,5]. Based on the data availability of vegetation growth and soil water content datasets, NDVI and GPP data were analyzed for 1982-2014, and VOD data for Table 2). GIMMIS-NDVI3g was used in the main text. Vegetation types were defined following the International Geosphere-Biosphere Program based on Terra and Aqua combined Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type (MCD12Q1) Version 6 data [6] and they were also used to exclude non-vegetation areas.

Drought indices
Vapor pressure deficit (VPD) and soil water content (SWC) are the two common indices of plant water availability, reflecting atmospheric and soil dryness, respectively [7]. Four climatic data sets were used to investigate the dynamics of VPD, including TerraClimate [8], the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis-Interim (ERA-Interim) [9], the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) [10], and the ECMWF Reanalysis Version 5 for land applications (ERA5-Land) [11] (Supplementary   Table 2). Two observation-driven soil moisture data sets were used to quantify the spatiotemporal dynamics of SWC. Since vegetation mainly absorbs water through the root system, the root zone soil moisture data as a proxy for total soil water content (SWC), was also used [12], which includes satellite-based soil moisture from GLEAMv3.5a [13][14][15] and the Global Land Data Assimilation System Version 2 (GLDAS-2) from the Noah Model 3.6 [16] (Supplementary Table   2). The two soil moisture products were generated from different data assimilation systems and different climatic forcing inputs. ERA-Interim VPD and GLEAM SWC were used.in the main text. All VPD and SWC data sets were resampled to a 1° spatial resolution. To ensure comparability among different observation-driven data sets, we transformed the average total soil moisture of the study domain ( Supplementary Fig. 1) as a fraction (%) of the simulated climatological mean during the baseline period 1961-1990 [17].

Aridity index
The aridity index (AI) measures the balance between the atmospheric water supply to the land (precipitation) and its demand from the land surface (potential evaporation, Ep). To consolidate the study findings and avoid bias from data uses, we used five precipitation products and three reanalysis-based Ep products (Supplementary Table 2), to estimate the AI and delineate drylands extent. We derived AI estimates by combining all possible pairs of precipitation and Ep data-based products. For the longer 1958-2014 period, we derived nine ensemble members of AI using three precipitation products (the Climatic Research Unit TS4.05 (CRUTS4.05) [18], the Princeton Global Meteorological Forcing v3 (PGFv3) [19], and the University of Delaware v5.01 (Udelv5.01) [20]) and three Ep products (GLDAS-2 [16], the National Center for Environmental Prediction/the National Center for Atmospheric Research (NCEP/NCAR) [21] and PGFv3 [19]).
For the shorter 1982-2014 period, we derived six members of AI using all available precipitation and Ep databases. All AI data sets were aggregated to a 1° spatial resolution. The reliabilities of the 15 members were verified by comparing each annual mean AI ( Supplementary Fig. 2).

Dynamic global vegetation models
We used 18 GPP and 14 leaf area index (LAI) simulations of offline dynamic global vegetation models (DGVMs) from the "Trends in net land atmosphere carbon exchange" (TRENDY-v9) project for 1982-2014 (Supplementary Table 3). For TRENDY DGVMs, we used S3 simulations that incorporated changes in climate forcing, rising atmospheric CO 2 concentrations, and land-use change. GPP and LAI for the vegetated land were resampled to a 1° spatial resolution, and aggregated to the growing season period (April to October) on a yearly timescale.

Earth system model outputs
We used monthly outputs of GPP, near-surface air temperature (T air ), near-surface relative humidity (RH), and total soil water content (SWC) from eleven Earth system models (ESMs) (Supplementary Table 4 were resampled to a 1° spatial resolution, and aggregated to the growing season period (April to October) on a yearly timescale.

Flux and environmental measurements
The carbon flux monitoring was carried out at Naqv alpine meadow ecosystem (AME) station (31.64°N, 92.01°E, 4600 m above sea level), which is a standard observation site of the China Flux Observation and Research Network. The station has a short growing season due to low temperatures, low SWC, and concentrated precipitation in Summer [22]. At this study site, Kobresia pygmaea accounts for about 70% of the total community aboveground biomass and 90% of the total community coverage. Thus, we used the phenological period of Kobresia Pygmaea to represent the growing season length (GSL) in AME. We randomly selected ten clusters of Kobresia Pygmaea in a circular area with a radius of 200 m around the flux tower, marked them with flags and observed them daily. We recorded the greening time of each cluster and calculated the average value as DOY at the start of the growing season. We recorded the total yellowing time of each individual plant and calculated the average value as DOY at the end of the growing season. For more details about flux data processing and phenological observation, see the Ref. [22].

Drylands extent with different aridity metrics
As a commonly used land surface aridity index, AI possesses high uncertainty in fully representing the complexity of aridity changes [17]. Meteorological conditions shape the geographic patterns of VPD and SWC at the land surface. We characterized how these aridity metrics had changed in recent decades (1958-2014; 1982-2014) as follows.
We first established an empirical relationship between a given aridity metric (VPD and SWC) and AI using all grid points in the AI-defined drylands (corresponding to AI = 0.65; Supplementary Fig. 1) [17,[23][24][25]. Considering potential non-linearity in the relationship, we tested three different empirical models: linear, polynomial, and exponential. The model with the maximum R 2 was adopted to describe the relationship between a given variable and AI [17]. With the pre-established regression models, we then calculated the respective thresholds of VPD and SWC for delineating dryland extent ( Supplementary Fig. 4), based on which year-to-year variations of drylands were displayed. The baseline period was set as 1961-1990 (or a subset of years thereof depending on data availability) [17]. The fraction indices of f atm and f soil , for VPD and SWC, respectively, represent percentage changes relative to the baseline, where a positive (negative) value implies an expansion (contraction) of drylands compared to the baseline period of 1961-1990.

Percentile binning
To investigate yearly vegetation response to water stress, the binning method was applied [7]. This method is capable of decoupling VPD and SWC changes, describing potential non-linearity, and observing correlated and causal relationships. For each pixel, we determined segmentation threshold values of the 20th, 40th, 60th, 80th, and 100th percentile of VPD and SWC, which will later be used to bin the data. Data of all variables (VPD, SWC, and vegetation growth indices) are sorted into 5 bins according to the 0-20th, 20-40th, 40-60th, 60-80th, and 80-100th percentiles of VPD or SWC. Next, within each SWC bin ( f = 1, 2, 3, 4, 5), the ranking from minimum ,min For example, VPD stress on NDVI after excluding the effects of VPD-SWC coupling (termed ΔNDVI (VPD|SWC)) was derived from the changes in NDVI from low VPD to high VPD at each SWC bin. We calculated the difference between the highest VPD and the lowest VPD corresponding to NDVI in each SWC bin (ΔNDVI (VPD|SWC)): where G is the total number of VPD bins, g is the specific VPD bin number, ,min g m and ,max g m is the minimum and maximum SWC bin number at the VPD bin g .
In the Liu et al.' (2020) [7] paper, the percentile bins were conducted using daily data because VPD and soil moisture are diagnosed to be decoupled at the daily or finer scales [7]. To further test the long-term response of vegetation to water stress, in addition to the logical function based on prior knowledge, we also examined the Spearman partial correlation between vegetation indices and SWC or VPD after statistically controlling changes in the other drivers. For example, the partial correlation coefficient between NDVI and SWC was calculated by excluding the influence of VPD. consequently saturated water vapor concentration over land (es, red bars at lower right) exceed growth in actual water vapor concentration (ea, blue bars). Increases in sensible and latent heat (associated, respectively, with temperature and water vapor, and represented by the area of each bar) have the same amount over land and ocean, with sensible heat increasing more over land than oceans and latent heat increasing more over oceans. Relative humidity (ratio of blue to red bar length) decreases over land [26]. Td: dew point temperature; Red cycle: air mass; Black arrow: moving direction. In more detail, atmospheric dryness is generally regarded as a simple thermodynamic consequence of warming [27]. The greater warming occurs over drylands than the ocean [27]. es increases more over drylands than over oceans because es is driven only by temperature. It impedes the transport of moist air masses from oceans to drylands. Oceanic air masses adverted over drylands contain insufficient water vapor to keep pace with the greater increase in es over drylands [26,28]. Hence, the increase of near-surface specific humidity over drylands is relatively low and contains insufficient ea to follow Clausius-Clapeyron scaling (~7% per ℃). The enlarged contrast between saturated and actual water vapor induced higher vapor pressure deficit (VPD). Figure 29. A simplified scheme linking soil water content (SWC) deficit to plant water stress.

Supplementary
The roles of climate, soil, and vegetation are linked to plant response through two fundamental processes: the first one is that SWC dynamics control the intensity and duration of SWC deficit periods, while the second one is that SWC deficit regulates plant physiology through the plant water potential, which in turn affects cell turgor and the relative water content of the plant's living cells under pressure [29].