Response of Vegetation to Drought in the Source Region of the Yangtze and Yellow Rivers Based on Causal Analysis

Abstract

The vegetation and ecosystem in the source region of the Yangtze River and the Yellow River (SRYY) are fragile. Affected by climate change, extreme droughts are frequent and permafrost degradation is serious in this area. It is very important to quantify the drought–vegetation interaction in this area under the influence of climate–permafrost coupling. In this study, based on the saturated vapor pressure deficit (VPD) and soil moisture (SM) that characterize atmospheric and soil drought, as well as the Normalized Differential Vegetation Index (NDVI) and solar-induced fluorescence (SIF) that characterize vegetation greenness and function, the evolution of regional vegetation productivity and drought were systematically identified. On this basis, the technical advantages of the causal discovery algorithm Peter–Clark Momentary Conditional Independence (PCMCI) were applied to distinguish the response of vegetation to VPD and SM. Furthermore, this study delves into the response mechanisms of NDVI and SIF to atmospheric and soil drought, considering different vegetation types and permafrost degradation areas. The findings indicated that low SM and high VPD were the limiting factors for vegetation growth. The positive and negative causal effects of VPD on NDVI accounted for 47.88% and 52.12% of the total area, respectively. Shrubs were the most sensitive to SM, and the response speed of grassland to SM was faster than that of forest land. The impact of SM on vegetation in the SRYY was stronger than that of VPD, and the effect in the frozen soil degradation area was more obvious. The average causal effects of NDVI and SIF on SM in the frozen soil degradation area were 0.21 and 0.41, respectively, which were twice as high as those in the whole area, and SM dominated NDVI (SIF) changes in 62.87% (76.60%) of the frozen soil degradation area. The research results can provide important scientific basis and theoretical support for the scientific assessment and adaptation of permafrost, vegetation, and climate change in the source area and provide reference for ecological protection in permafrost regions.

1. Introduction

With the joint impact of climate change and intense human activities, complex extreme drought events are becoming more and more frequent, and the problem of ecosystem imbalance under water constraints is more acute [1]. The source region of the Yangtze River and the Yellow River (SRYY) is an important ecological security barrier in China and even Asia, which is sensitive to climate change. The average warming rate in the source region from 1961 to 2021 was 0.36 °C/10a. Although the region as a whole is humid and warm, local extreme droughts still occur frequently [2]. Water status is the biggest limiting factor affecting vegetation growth. Long-term continuous drought events affect the physiological state of vegetation in the source area through changes in factors such as temperature, precipitation, and soil water content, thereby affecting the carbon cycle and vegetation greenness [3]. Among them, vapor pressure deficit (VPD) and soil moisture (SM) have received extensive attention in exploring the mechanism of water constraints on vegetation growth. On the one hand, plant death may be due to the imbalance of hydraulic transport process between ‘soil-plant system’ caused by xylem catheter embolism; on the other hand, high VPD-induced stomatal closure leads to decreased photosynthesis, which makes the photosynthetic products insufficient to support the normal physiological metabolic processes of plants [4]. In addition, the source area is widely distributed with a large area of permafrost and seasonal frozen soil, and the surface experiences a freeze–thaw process for more than half a year. The freeze–thaw cycle changed the regional water cycle and the characteristics of recharge, runoff and discharge, resulting in a profound change in the interconnection between drought and alpine vegetation in the SRYY [5]. Due to the different biological characteristics and ecological adaptability of different vegetation types, the response to drought is also different (Figure 1). Under the influence of climate–frozen soil coupling, clarifying the mechanism of VPD and SM on vegetation continues to be a challenging endeavor.

At present, most research focuses on analyzing the change characteristics of drought and vegetation or exploring the correlation [6,7]. For drought monitoring and simulation, according to the different causes, drought can be divided into atmospheric drought, soil drought and physiological drought. Atmospheric drought caused by insufficient precipitation will have different degrees of impact on soil, vegetation, runoff and groundwater, which in turn will cause agricultural drought and ecosystem drought [8]. Drought stress can cause water deficit and osmotic stress in plants, which will affect vegetation phenology and growth and development, and even lead to plant death in severe cases [9]. Many scholars have carried out research on the response of vegetation to drought, which has revealed the significant impact of drought on vegetation productivity. In addition, the response of vegetation to drought is very complex, which is closely related to vegetation type [10,11], dry and wet environment and drought time scale [12].

Traditional remote sensing vegetation indices such as Normalized Differential Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) have been widely used in large-scale vegetation drought monitoring to analyze the impact of drought on vegetation [13]. However, the above indexes represent the greenness and canopy structure of vegetation and lack direct connection with vegetation photosynthesis. When vegetation is subjected to water stress, they cannot reflect the short-term changes of vegetation photosynthesis caused by drought stress in time. A field experiment showed that when the vegetation was subjected to short-term severe drought stress, NDVI remained at a high level, and there was a significant lag phenomenon [14]. Selecting sensitive indicators to detect the response of vegetation to drought is helpful to grasp the situation of vegetation affected by drought as early as possible, and to help understand the response process of vegetation to drought stress. A large number of studies have shown that SIF can quickly reflect the effect of water stress on the physiological growth of vegetation [15]. SIF is the light signal of red and far-red bands released by chlorophyll in the excited state after capturing light quantum in the process of photosynthesis, which is directly related to the photosynthetic rate of vegetation. It is an ‘associated product’ of vegetation photosynthesis. With the advancement of remote sensing technology, there are a series of mature SIF products [16]. GOSIF is a long-term sequence data with higher spatial resolution and global continuous coverage developed by integrating discrete SIF detected by OCO-2, MODIS remote sensing data and meteorological reanalysis data [17]. It has the advantages of long time series and high spatial resolution and can be used to study the response of vegetation to drought.

In addition, it was pointed out that high VPD and low SM were two types of water stress in vegetation. However, which type of drought plays a leading role in vegetation growth is a highly controversial topic [18]. Several studies have indicated that the increasing atmospheric drought since 1990 has greatly weakened the carbon absorption capacity of global ecosystems, offsetting some of the carbon sequestration caused by the CO₂ fertilization effect [19]. However, Liu et al. [20] used the Binning method to decouple SM and VPD and indicated that in over 70% of the global regions, SM is a major determinant of the growth of vegetation. This is particularly true for semi-arid ecosystems. It can be seen that there are differences in the research results of vegetation response to drought, and the mutual feedback mechanism is still unclear. The coupling relationship between VPD and SM has obvious spatial heterogeneity. Studies have shown that soil–atmosphere composite droughts are projected to be more frequent and extreme [21], which will strongly limit the carbon absorption capacity of vegetation. Therefore, it is of great significance to clarify the impacts of VPD and SM on evolution of vegetation. In addition, the identification of the interaction between the two involves the causal relationship between hydrological and ecological factors. The related research is changing from correlation analysis to causal analysis. The causal analysis algorithm [22] and structural equation model [23] have also been gradually applied to quantify the impact of environmental driving factors on evolution of vegetation. Therefore, this study is supported by causality discovery technology to identify the causal interaction information between vegetation and drought.

The plateau ecosystem in the SRYY is fragile, and the grassland degradation in some areas affected by drought events is severe [24]. The coupling effect of atmospheric–soil drought in the source area affected by permafrost landform and high-altitude climate has regional particularity, and the spatial heterogeneity of the ecosystem is obvious [25]. In view of this, the main aims of this study include (1) revealing the temporal and spatial evolution and variation trend of land–atmosphere water constraint events and vegetation greenness and productivity in the SRYY; (2) the causal analysis chain between Normalized Differential Vegetation Index (NDVI), Solar-induced fluorescence (SIF) and VPD; SM is constructed to explore the mechanism of the land–atmosphere water constraint on vegetation growth; (3) zoning (permafrost degradation zone) classification (different vegetation types) to analyze the interaction of vegetation–drought.

2. Materials and Methods

2.1. Study Area

The SRYY is situated in the eastern margin of the Qinghai-Tibet Plateau, between 90.55°E~103.41°E, 32.16°N~36.12°N. The water system in the area is developed, the river network is dense, and there are many glaciers and lakes. It is the birthplace and catchment area of the Yangtze River and Yellow River [26]. The area of the source area is about 26.44 × 104 km², and the terrain environment is complex and diverse, with an average altitude of 4464 m (Figure 2a,b). Affected by the temperature rise, although the drought degree in the source area has alleviated, extreme droughts still occur frequently. For example, the most severe drought occurred in 2017, reaching moderate to severe drought, and 243.9 khm² of crops were affected. Affected by climate change, the transformation relationship between water cycle elements in the region will be more complex in the future, and the complexity and uncertainty of drought risk will be further aggravated [27].

Due to the control of altitude climate and permafrost landform, the composition of vegetation flora in the source area is simple. The predominant vegetation in the source area is grassland, encompassing approximately 76% of the entire area. In addition, there are scattered forest ecosystems composed of forest land, sparse forest land and a small amount of shrub land (Figure 2c). There are various types of frozen soil in this area. In this study, the surface freezing number model was used to divide the boundary of the permafrost degradation area, constituting 18.71% of the total source area and predominantly focused in the southern region of the source region of the Yangtze River (SRYA) (Figure 2d).

2.2. Data

The data consist of meteorological data, soil moisture data, vegetation data, vegetation type data, etc. The data information and its sources are detailed in

Table 1. The information of data.

Type	Data	Information	Source
meteorological data	2 m temperature	Monthly data from 1980 to 2021, data format is nc and spatial resolution is 0.1° × 0.1°.	ERA5 (https://cds.climate.copernicus.eu, accessed on 15 September 2023)
	2 m dewpoint temperature	Monthly data from 1980 to 2021, data format is nc, and the spatial resolution is 0.1° × 0.1°.
	temperature	Hourly data from 2000 to 2020, data format is nc, and the spatial resolution is 0.1° × 0.1°.
soil data	SM	Monthly data from 1980 to 2021, data format is nc, and spatial resolution is 0.5° × 0.5°.	GLEAM (https://www.gleam.eu/, accessed on 15 September 2023)
vegetation data	solar-induced fluorescence	Monthly data from 2001 to 2021, spatial resolution is 0.05° × 0.05°.	Global Ecology Group Data Repository: (http://globalecology.unh.edu/data/GOSIF.html, accessed on 15 September 2023)
	Normalized Differential Vegetation Index	16 d data from 2001 to 2021, spatial resolution is 250 m × 250 m.	NASA (https://modis.gsfc.nasa.gov/, accessed on 15 September 2023)
vegetation type data	GLC_FCS30-2020	Global 30 m fine land cover	https://data.casearth.cn, accessed on 26 September 2023

.

2.2.1. Meteorological Data

ERA5 (ECMWF Reanalysis V5) is the fifth-generation global reanalysis product produced by the European Centre for Medium-Range Forecasts (ECMWF). This study used monthly and hourly data from 1980 to 2021 with a spatial resolution of 0.1° × 0.1°, including 2 m air temperature and 2 m dew point temperature (https://cds.climate.copernicus.eu/, accessed on 15 September 2023). The data set is mainly used to calculate VPD. After calculating the VPD, in order to cooperate with the subsequent calculation and unify the resolution, it is resampled to 250 m × 250 m.

2.2.2. Soil Data

The Global Land Evaporation Amsterdam Model (GLEAM) is a set of algorithms for the estimation of terrestrial evaporation and root-zone SM from satellite data (www.GLEAM.eu, accessed on 15 September 2023). These output datasets can be utilized for hydrological applications, climate studies and research about land–atmosphere feedbacks. In this study, we used global SM observation data of each month from 1980 to 2021 and interpolated the spatial resolution to 250 m × 250 m by using the bi-linearly interpolated method.

2.2.3. Vegetation Data

In this study, the NDVI data were obtained from the MOD13A3 dataset, which was supplied by the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) National Aeronautics and Space Administration (NASA) EOSDIS Land Processes Distributed Active Archive Center (https://appeears.earthdatacloud.nasa.gov/, accessed on 15 September 2023). The data consist of 483 images from 2001 to 2021, with a spatial resolution of 250 m and a temporal resolution of 16 days. A study showed that the one-year NDVI calculated by the Maximum Value Composite (MVC) method could reflect the potential productivity of vegetation. In contrast, the one calculated by the arithmetic mean (AM) method can reflect vegetation’s annual average production status. Therefore, we used the AM method to process 23 images in each natural year, generating the annual mean values of the NDVI. It will be used to analyze the spatiotemporal variations of vegetation and its response to climatic conditions.

Li et al. [28] developed a high spatial resolution and spatial continuous long time series GOSIF data based on OCO-2 SIF data, MODIS remote sensing data and meteorological reanalysis data. The resolution of the data is 0.05° × 0.05°. Using the 2001–2021 GOSIF monthly scale product (http://globalecology.unh.edu/, accessed on 15 September 2023), in order to match the spatial resolution of other data, the causal analysis calculation is resampled to 250 m × 250 m.

2.2.4. Land Cover Data

The global land-cover product with fine classification system (GLC_FCS) is a product containing 30 land cover types generated by Zhang et al. [29] using multi-temporal Landsat spectral and texture features and corresponding training data. Land cover types include rainfed/irrigated farmland, evergreen/deciduous broad-leaved forest, evergreen/deciduous coniferous forest, mixed forest land, shrub, grassland, wetland, artificial surface, tundra, sparse vegetation, bare land, water body, permanent ice and snow. GLC_FCS is a global land cover product with a spatial resolution of 30 m.

2.3. Methods

According to the framework of ‘evolution characteristics-mutual feedback mechanism-partition statistics’, the key content and technical route of this research are put forward. Firstly, based on the VPD and SM indexes that characterize atmospheric drought and soil drought, the evolution of atmospheric and soil drought events in the SRYY was revealed. NDVI and SIF were selected to analyze the characteristics and trends of vegetation greenness and productivity in the source area. Secondly, by constructing the causal analysis chain between NDVI/SIF and VPD/SM, the mechanism of land–atmosphere water constraint on vegetation growth was explored, the dominant factors affecting vegetation change were identified, and the lag effect of vegetation on drought was quantified. In addition, the response of permafrost degradation areas and different vegetation types to water constraints was analyzed in depth. The technical route is shown in Figure 3.

2.3.1. Division Method of Permafrost Degradation Zone

The surface freezing number model uses the ratio of the surface freezing index to the thawing index to divide the boundaries of different frozen soil types [30].

$$FI=\sum _{i=1}^{NF}\left|T_i\right|,T_i<0 °\mathrm{C}$$

(1)

$$TI=\sum _{i=1}^{NT}\left|T_j\right|,T_j\ge 0 °\mathrm{C}$$

(2)

In the formula: FI, TI and T_i are the freezing index, the melting index, and the average temperature of the i and j months, respectively.

$$\mathrm{F}=\frac{\sqrt{FI}}{\sqrt{FI}+\sqrt{TI}}$$

(3)

0 ≤ F < 0.5 indicates that there is no permafrost, and F ≥ 0.5 indicates that there is permafrost.

2.3.2. Calculation of VPD

VPD is the difference between e_s and e_a [31]. The formula is as follows:

$$VPD=e_s-e_a$$

(4)

The saturated vapor pressure (e_s):

$$e_s=6.11\times \exp (17.27\times (T+237.3))$$

(5)

The actual vapor pressure (e_a):

$$e_a=6.11\times \exp (17.27\times (T_{dew}+237.3))$$

(6)

where T is 2 m temperature, T_dew is 2 m dewpoint temperature, and the unit is °C.

2.3.3. Sen–Mann–Kendall

Calculate the statistical value of N to the data, expressed as

$$\beta =\frac{(X_i-X_k)}{\left(i-k\right)}$$

(7)

In the formula, X_i and X_k represent the data values at time i and k, respectively. For time series, the standardized test statistic Z is defined. The method for judging the significance of the Sen–Mann–Kendall trend [32] is shown in

Table 2. Trend division standard.

β	Z	Type	Trend Features
β > 0	2.58 < Z	4	Extremely significant increase
	1.96 < Z ≤ 2.58	3	Significantly increased
	1.65 < Z ≤ 1.96	2	Micro-significant increase
	Z ≤ 1.65	1	Not significantly increased
β = 0	Z = 0	0	No change
β < 0	Z ≤ 1.65	−1	Not significantly reduced
	1.65 < Z ≤ 1.96	−2	Micro-significant reduction
	1.96 < Z ≤ 2.58	−3	Significantly reduced
	2.58 < Z	−4	Extremely significant decrease

.

2.3.4. Peter–Clark Momentary Conditional Independence (PCMCI)

The identification of the interaction between drought and vegetation involves the causal relationship of multiple factors, which needs the support of causal relationship discovery technology. This technique uses prior knowledge contained in time series measurement and variable selection to extract causal interaction information between variables in a given system. Different from the traditional correlation method, the correlation method has no way to determine the direction of causal connection and can only obtain an undirected graph. The purpose of the causal analysis method is to find the cause of a certain result so as to construct a directed causal network structure between system variables [33]. The transformation from correlation analysis to causal analysis is an inevitable trend of research and development.

The theories and algorithms of causal inference have made significant progress in the past few decades. Granger made important work in the early empirical causal inference methods [34]. Based on the theory of causality discovery, the algorithms of causality discovery theory are systematically sorted out from the aspects of basic assumptions, principles, advantages and disadvantages (

Table 3. Comparison of various causal inference algorithms.

Range	Fundamental Assumption	Principle	Advantages and Disadvantages
GC (Granger Causality)	The cause precedes the result; linear relationship between variables	vector autoregression	High-order condition sets appear at runtime, which reduces the effectiveness of the algorithm.
TE (Transfer Entropy)	The cause precedes the result; reasons provide useful information for the prediction of results.	conditional mutual information	It requires the probability distribution of the variables; same as above.
PC (Peter–Clark algorithm)	The cause precedes the result	graphics	Intuitive; does not require the determination of high-order independence
CCM (Convergent Cross Mapping)	Primitive manifold and reconstruction shadow manifold	/	It cannot be used for time series with strong coupling relationship.

). It is found that the PC algorithm is more intuitive than other algorithms, does not require the determination of high-order independent connections and has higher operating efficiency.

In view of the fact that the results of the original PC algorithm may contain false positive causality, in this study, the PC algorithm superimposed with the MCI test (PCMCI algorithm) was used to study the interaction. The PCMCI algorithm is an improvement of the PC algorithm. Based on the dependence on the lag time, the algorithm can efficiently process high-dimensional data and eliminate the false positive causality between variables so as to obtain a more credible causal relationship between variables. The PCMCI algorithm mainly includes two steps [35]:

(1) PC₁ condition selection: the PC algorithm is used to estimate the superset of the parent variable $\mathrm{P}({\mathrm{X}}_{\mathrm{t}}^{\mathrm{j}})$ of all variables, that is, $\widehat{\mathrm{P}}({\mathrm{X}}_{\mathrm{t}}^{\mathrm{j}})$;

(2) Take the parent variable $\mathrm{P}({\mathrm{X}}_{\mathrm{t}}^{\mathrm{j}})$ in step (1) as the condition. The momentary conditional independence test (MCI) algorithm is used to test the variable pairs ‘$({\mathrm{X}}_{\mathrm{t}-\mathsf{\tau }}^{\mathrm{i}},{\mathrm{X}}_{\mathrm{t}}^{\mathrm{j}})$’ (where ‘$\mathrm{i},\mathrm{j}\in \left\{1,2,\dots ,\mathrm{N}\right\}$’}, time delay ‘$\mathsf{\tau }\mathsf{ϵ}\left\{1,2,\dots ,{\mathsf{\tau }}_{\max }\right\}$’). And establish the connection from ${\mathrm{X}}_{\mathrm{t}-\mathsf{\tau }}^{\mathrm{i}}\to {\mathrm{X}}_{\mathrm{t}}^{\mathrm{j}}$, which is expressed as

$$\mathrm{M}\mathrm{C}\mathrm{I}={\mathrm{X}}_{\mathrm{t}-\mathsf{\tau }}^{\mathrm{i}}⫫{\mathrm{X}}_{\mathrm{t}}^{\mathrm{j}}\left|\widehat{\mathrm{P}}\left({\mathrm{X}}_{\mathrm{t}}^{\mathrm{j}}\right)\backslash \{{\mathrm{X}}_{\mathrm{t}-\mathsf{\tau }}^{\mathrm{i}}\right\},\widehat{\mathrm{P}}\left({\mathrm{X}}_{\mathrm{t}-\mathsf{\tau }}^{\mathrm{i}}\right)$$

(8)

The above two steps, condition selection and the MCI algorithm, both contain conditional independent tests, in which the test statistics involved use partial correlation test statistical methods. Based on the PCMCI algorithm, the causal influence value between the elements (represented by the maximum value of the absolute value of MCI) and its corresponding delay time can be obtained to form the relevant spatial grid data. The statistical values of causal effects and lag duration in different vegetation types and permafrost degradation areas were obtained by using regional statistical tools. According to the causal influence value between variables estimated by the PCMCI algorithm, the grading standard is shown in

Table 4. Causal impact value grading standard.

Level	Range	Causal Impact Strength	Level	Range	Causal Impact Strength
1	0.8 < x ≤ 1	Extremely strong positive causal impact	7	−0.2 < x < 0	Weak negative causal effect
2	0.6 < x ≤ 0.8	Relatively strong positive causal effect	8	−0.4 < x ≤ 0.2	Moderate negative causal effect
3	0.4 < x ≤ 0.6	Strong positive causal effect	9	−0.6 < x ≤ −0.4	Strong negative causal effects
4	0.2 < x ≤ 0.4	Moderate positive causal effect	10	−0.8 < x ≤ −0.6	Relatively strong negative causal effects
5	0 < x ≤ 0.2	Weak positive causal effect	11	−1 ≤ x ≤ −0.8	Extremely strong negative causal effects
6	x = 0	No causal effect

.

3. Results

3.1. Evolution of Drought and Vegetation in the SRYY

3.1.1. The Evolution Characteristics of VPD and SM

From 1980 to 2020, the annual average VPD in the source region of the Yangtze River and the Yellow River showed a fluctuating increasing trend. Compared with VPD, the variation of SM in the source region was weak, showing a trend of decreasing first and then increasing around 1991. In 1987, 2003 and 2021, the source area showed high VPD and low SM, while in 1983 and 2012, it was the opposite, indicating that VPD and SM were significantly negatively correlated. On the whole, the atmospheric drought in the source area has increased in the past 40 years, and the soil tends to be humid. Spatially, the spatial heterogeneity of atmospheric and soil drought is obvious. The VPD in the northeastern corner of the source region of the Yellow River (SRY) and the western part of the SRYA is higher than that in other regions, showing strong atmospheric drought. The SM values in the western part of the SRYA and the western part of the SRY are lower, showing obvious soil drought (Figure 4a–d).

Spatially, the VPD of 92.48% of the source area exhibited an upward trend. The degree of atmospheric drought in the eastern part of the SRY increased significantly, and the drought gradually weakened from east to west, and 69.16% of the region showed no significant increase (Figure 5a). On the contrary, the SM in 96.09% of the source area exhibited an upward trend, and more than half of the area had obvious soil wetting, which gradually increased from south to north (Figure 5b). In general, from 2000 to 2021, the degree of atmospheric drought in the SRYY increased, and the soil drought eased.

3.1.2. Evolution Characteristics of NDVI and SIF

The distribution of SIF and NDVI in the SRYY from 2001 to 2021 is illustrated in Figure 6a–d. The fluctuation ranges of annual average NDVI and SIF in the source area were 0.247–0.287 and 0.025–0.029, respectively. From the change trend of NDVI and SIF, the vegetation in the source area experienced an evolution process of ‘increase-decrease-increase’. Both of them indicated that the vegetation coverage improved significantly after 2008, and the vegetation improvement was the most obvious in 2010 and 2020. In the space, the vegetation change in the SRYY exhibits significant regional differences, with high performance in the southeast and low performance in the northwest. The spatial distribution of NDVI and SIF is very consistent, indicating that SIF can better reflect the true state of vegetation growth.

From 2001 to 2021, NDVI and SIF in the region exhibited an overall upward trend, and the area of vegetation improvement accounted for 82.26% and 63% of the total area, respectively (Figure 7). With the establishment of the nature reserve in 2001, the growth rate and area of vegetation have increased significantly, but there is still a downward trend in some areas. The distribution of SIF and NDVI changes has a good consistency. The areas with improved vegetation are mainly distributed in the northern part of the SRYA and the SRY, while the degraded areas are mainly distributed in the eastern part of the SRYA and the southern part of the SRY with more human activities.

3.2. Interaction Mechanism between Drought and Vegetation Based on Causal Analysis

3.2.1. Effects of Drought and Vegetation on the SRYY

The average causal effect of VPD on SM in the SRYY is −0.5426, and the overall performance of the two is a strong negative causal relationship. In the space, about 91.75% of the source regions are significantly correlated, of which 87.05% are negatively correlated with SM, and only a few regions in the SRYA are positively correlated, indicating that the coupling between VPD and SM in the source region is strong. High VPD is usually accompanied by low SM events. The average causal effect of VPD on NDVI in the source area was −0.022, and the causal effect of SM on NDVI was 0.1068. According to the classification standard of causal influence value, the positive and negative causal effects of VPD on NDVI accounted for 47.88% and 52.12% of the total area, respectively, with moderate and strong causal effects. SM was dominated by moderate or strong positive causal effects, accounting for 24.36% and 23.46% of the total area, respectively. It was mainly distributed in the west of the SRYA and on both sides of the SRY, showing that the vegetation became green with the increase in soil moisture (Figure 8a,

Table 5. The proportion of different causal effects between drought and NDVI in the source area (%).

Type	Extremely Strong (−)	Relatively Strong (−)	Strong (−)	Moderate (−)	Weak (−)	Weak (+)	Moderate (+)	Strong (+)	Relatively Strong (+)	Extremely Strong (+)
VPD-SM	4.03	47.56	35.47	7.26	0.20	0.04	0.76	2.42	2.21	0.08
SM-NDVI	0.13	3.68	12.49	17.32	4.68	5.34	24.36	23.46	8.14	0.4
VPD-NDVI	0.05	3.42	16.27	25.32	7.06	6.82	23.96	14.03	2.98	0.09
SM-VPD	4.02	44.18	32.83	7.07	0.12	0.09	1.87	7.17	2.65	0.01
NDVI-VPD	0.06	4.15	23.33	30.97	6.61	4.80	17.74	10.28	1.98	0.06
NDVI-SM	0.12	3.44	10.35	11.76	2.81	4.28	25.90	30.03	10.84	0.48

).

Soil and atmosphere participate in the water cycle and energy cycle on the regional scale, and have complex coupling feedback with vegetation. Through causal analysis, it was found that the causal effect of SM on VPD in the source area was −0.4764; that is, when the soil moisture decreased, the evapotranspiration decreased, resulting in a decrease in water vapor content in the atmosphere, a corresponding increase in VPD, and an increase in atmospheric drought. In the space, 65.13% of the source area is negatively correlated, mainly with medium negative causality, and the positive correlation area is mainly distributed in the middle of the SRYY. The average causal effects of NDVI on VPD and SM were −0.1257 and 0.1935, respectively. According to the classification standard of causal influence value, the area of negative causal influence of NDVI on VPD is 88.22%, which is mainly strong and strong negative. The strong positive causal effect on SM is dominant. In most areas, the better vegetation coverage is more conducive to increasing soil moisture and alleviating atmospheric drought (Figure 8b, Table 5).

For the impact of drought on vegetation in the source area, the average causal effect of VPD on SIF was −0.0775, and the causal effect of SM on SIF was 0.2803. According to the classification standard of causal influence value, the positive and negative causal effects of VPD on SIF accounted for 45.89% and 54.11% of the total area, respectively. And the influence area of atmospheric drought on SIF is consistent with NDVI. The strong positive causal effect of SM on SIF accounted for 31.70% of the total area. Compared with the influence of drought on NDVI, the influence of drought on SIF is slightly higher and the influence range is wider (Figure 9a). For the feedback of vegetation to drought, the average causal effects of SIF on VPD and SM were −0.0968 and 0.2359, respectively. According to the classification standard of causal influence value, the positive and negative causal effects of SIF on VPD accounted for 42.56% and 57.44% of the total area, respectively. The strong positive causal effect on SM is dominant. Compared with the feedback of NDVI to drought, the feedback of SIF to atmospheric drought is weaker, and the feedback to soil moisture is stronger. In summary, the effect of SM on vegetation greenness and productivity is stronger than that of air drought (Figure 9b,

Table 6. The proportion of different causal effects between drought and SIF in the source area (%).

Type	Extremely Strong (−)	Relatively Strong (−)	Strong (−)	Moderate (−)	Weak (−)	Weak (+)	Moderate (+)	Strong (+)	Relatively Strong (+)	Extremely Strong (+)
VPD-SM	4.02	46.51	33.18	8.32	0.24	0.04	1.51	3.19	2.83	0.14
SM-SIF	0.004	1.30	8.13	12.45	2.79	2.82	19.36	31.71	20.35	1.08
VPD-SIF	0.80	10.95	22.11	17.05	3.20	3.65	22.59	16.57	3.02	0.05
SM-VPD	3.98	43.31	30.31	8.13	0.08	0.04	3.19	8.22	2.71	0.04
SIF-VPD	0.07	6.36	24.40	20.92	5.73	7.18	23.11	10.75	1.51	0.03
SIF-SM	0.003	0.84	5.76	11.84	4.71	6.10	31.98	30.06	8.59	0.11

).

3.2.2. The Lag Time between Drought and Vegetation in the SRYY

The average delay time of VPD’s impact on SM in the SRYY is 0.22 months. In the space, the proportion of regions with the delay time of 0, 1 and 2 months for the causal effect of source VPD on SM was 82.14%, 11.87% and 5.99%, respectively. Soil moisture in most regions was rapidly affected by atmospheric drought, and only the central and southwestern source areas were delayed by more than one month. The average delay time of the influence of VPD and SM on NDVI in the source area was 0.8479 and 0.7529, respectively, and the impact of SM on vegetation was rapid. In the space, the proportion of regions with the delay time of 0, January and February for the causal effect of SM on NDVI was 47.43%, 29.85% and 22.72%, respectively. The proportion of regions with delay time of 0, 1 and 2 months of VPD on NDVI was 40.74%, 33.72% and 25.54%, respectively. The impact of atmospheric drought on vegetation in the SRY was more rapid than that in the SRYA (Figure 10a).

The average delay time of the effect of source SM on VPD is 0.2285, and the feedback of SM on VPD is slightly slower than the effect. In the space, the proportion of regions with the delay time of 0, 1 and 2 months for the causal effect of source SM on VPD was 87.01%, 3.61% and 9.37%, respectively. The average delay time of VPD and SM in the source area to NDVI feedback was 0.9171 and 0.9377, respectively, and vegetation had a rapid effect on soil moisture. In the space, the proportion of regions with the delay time of 0, 1 and 2 months for the causal effect of NDVI on SM was 42.17%, 23.96% and 33.88%, respectively. The proportion of regions with the delay time of 0, 1, and 2 months for the causal effect of NDVI on VPD was 38.27%, 30.10% and 31.63%, respectively (Figure 10b).

For the impact of drought on vegetation in the source area, the average delay time of VPD and SM on SIF was 0.9233 and 0.9435, respectively, and the impact of atmospheric drought on SIF was relatively rapid. In the space, the proportion of regions with a delay of 0, 1 and 2 months for the causal effect of SM on SIF was 20.01%, 64.46% and 15.52%, respectively. The regions with a delay of 1 month were mainly distributed in the SRYA. The proportion of regions with a delay time of 0, 1 and 2 months for the causal effect of VPD on SIF was 27.17%, 52.17% and 20.66%, respectively. The regions with a delay of 1 month were mainly distributed in the SRY. The effect of drought on NDVI was faster than that of SIF (Figure 11a). For the feedback of vegetation to drought, the average delay time of VPD and SM to SIF feedback was 0.8437 and 1.0196, respectively, and vegetation had a rapid impact on atmospheric drought. In the space, the proportion of regions with 0, 1, and 2 months of delay in the causal effect of SIF on SM and VPD is roughly the same. Compared with the feedback of NDVI to drought, the feedback of SIF to atmospheric drought is faster, and the feedback to soil moisture is slower. In summary, except for the feedback of VPD to SIF, the interaction between drought and NDVI was faster than that of SIF (Figure 11b).

3.3. Response Mechanism of Vegetation to Drought

3.3.1. Responses of Different Vegetation Types to Drought

There are notable variations in the sensitivity of different vegetation types to drought response. For soil drought, the basic performance is a positive effect. The MCI coefficients of shrubs and NDVI/SIF were the largest, reaching 0.24 and 0.48, respectively, followed by sparse vegetation, grassland and cultivated plants, which had little effect on coniferous forests and broad-leaved forests. This shows that shrubs and herbs are sensitive to soil drought, and forest vegetation has a strong ability to resist drought stress. For atmospheric drought, drought–vegetation basically showed a negative effect. Cultivated plants were the most sensitive to atmospheric drought, and the MCI values of NDVI/SIF and VPD were −0.20 and −0.21, respectively, followed by shrubs, coniferous forests, broad-leaved forests, sparse vegetation and grassland. Compared with NDVI, all kinds of vegetation SIF have more obvious response to atmospheric drought. It shows that the photosynthesis of forest vegetation is more sensitive to atmospheric drought than that of herbs (Figure 12a). The response of NDVI to SM from fast to slow was cultivated plant > coniferous forest > grassland > broad-leaved forest > sparse vegetation > shrub; that of SIF to SM was cultivated plants > shrubs > sparse vegetation > grassland > coniferous forest > broad-leaved forest. The basic performance is that the forest has a long lag period for soil drought. The lag effect of vegetation on VPD is opposite, the forest response is faster, the lag period of herbaceous vegetation is longer, and the forest NDVI has a mainly positive response (Figure 12b).

3.3.2. Response Mechanism of Vegetation to Drought in Permafrost Degradation Area

The increase in active layer depth caused by permafrost degradation will lead to the redistribution of water and heat, thus affecting soil moisture and vegetation growth. Therefore, the response degree and lag time of vegetation to drought in the permafrost degradation area are extracted, as shown in Figure 13. The average causal effects of NDVI and SIF on SM were 0.21 and 0.41, respectively. In the space, 71.8% and 85.22% of the regions were positively correlated, and the whole was dominated by strong/relatively strong positive causality. The average causal influence values of VPD with NDVI and SIF were only −0.07 and −0.03. Therefore, due to the influence of permafrost degradation, the effect of SM on vegetation is more obvious. In addition, the response lag time of NDVI and SIF to SM in the permafrost degradation area was 23 days and 29 days, respectively. In the space, the lag time of NDVI and SIF in most areas was 1–2 months. The lag time of NDVI and SIF to VPD was 28 days and 31 days, respectively. Therefore, the response of vegetation to SM in the permafrost degradation area is faster and more sensitive.

Figure 14 shows the main factors that have causal effects on NDVI and SIF. On the whole, 60.63% of the source region SM has a stronger causal impact on NDVI than VPD, and 39.37% of the source region VPD has a stronger causal effect on NDVI than SM. Both are evenly distributed in space. The causal impact of SM on SIF was stronger than that of VPD in 60.04% of the source region, and the causal impact of VPD on SIF was stronger than that of SM in 39.96% of the source region. Spatially, the areas dominated by VPD are mainly concentrated in the SRY, and the areas dominated by SM are mainly distributed in the SRYA. For the permafrost degradation area, SM dominated NDVI (SIF) changes in 62.87% (76.60%) of the area. In summary, regardless of the source area or the permafrost degradation area, SM plays a significant role in the growth of vegetation in a wider range, and the permafrost degradation area is more obvious.

4. Discussion

4.1. Drought Evolution Law in the SRYY

As the ‘sensitive area’ and ‘start-up area’ of climate change, the uncertainty of extreme hydrological events in the SRYY has increased and has a wide and profound impact on terrestrial ecosystems. Soil and atmospheric drought are two categories of water stress experienced by vegetation. This study discovered an upward trend in the VPD of the source region, aligning with earlier findings [36]. VPD is controlled by multiple variables such as surface temperature and water vapor pressure, and the rapid temperature rise induces the increase of atmospheric saturated vapor pressure. However, the actual water vapor pressure did not rise at a consistent rate, forcing the difference between the two to gradually expand [37]. He et al. [38] also pointed out that with the intensification of climate change in the future, atmospheric drought may lead to systemic risks of climate. The SM in the source area showed an increasing trend, and the spatial distribution heterogeneity was strong, which may be related to the complex terrain and vegetation distribution in the source area. The vegetation coverage in the central and western parts of the SRYY is low, mainly grassland, and the storage capacity of soil water is insufficient, resulting in a low level of SM.

Studies have shown that soil–atmosphere composite droughts are projected to be more frequent and extreme [39]. The coupling effect of atmosphere–soil drought in the source area is strong, and 99.38% of the area is negatively correlated. Granger causality test showed that VPD and SM showed a two-way causal relationship in most areas of grassland in China [40]. Low SM will reduce land evapotranspiration, resulting in higher temperature and VPD; elevated VPD promotes the enhancement of land evapotranspiration, which in turn forces the dissipation of soil moisture [41]. However, after the frequency statistics of composite drought events in the source area, it is found that the composite drought in the source area shows a slowing trend from 1980 to 2020. The probability of composite drought in most areas of the source area is low, and only 33% of the total area of the SRYY occurs more than 10 times. Different from other regions, due to the unique geographical and geological environment and changeable climatic conditions in the SRYY, the soil moisture is easily affected by the complex multiphase water cycle such as the melting of glaciers, snow and frozen soil. The composite drought in the SRYY is concentrated in the northwest of the SRYA. The composite drought has a more profound and serious impact on the ecosystem than the single water-constraint event. In the future, the risk prevention and emergency response of the composite drought disaster in this region should be strengthened (Figure 15).

4.2. Comparison of NDVI and SIF Performance Differences in SRYY

This study shows that in the past 21 years, 82.26% of the regional vegetation greenness in the source area has increased significantly, but there is still a deterioration trend in some areas, and the vegetation degradation in the SRYA is more obvious. It is consistent with the research conclusion of Zhang et al. [42]. The research of Myers-Smith et al. [43] and Choler et al. [44] shows that with global warming, the vegetation index in the Arctic and Alpine regions is on the rise as a whole, indicating that climate change has a positive effect on the alpine vegetation. Climate change is the decisive factor affecting vegetation growth, but human protection and management actions can promote the rate of vegetation improvement in the short term. This study found that the vegetation in the source area improved significantly after 2008, and the overall vegetation coverage was characterized by high southeast and low northwest. On the one hand, the human disturbance in the source area also showed an increasing trend from northwest to southeast, the grassland area was under the intensity of severe human disturbance on average, and the spatial distribution of the two was similar [45]. The southeast is the key area for the implementation of ecological projects. The intensity of human interference is large and the vegetation coverage is high, which verifies the positive impact of the execution of ecological restoration in the southeast. On the other hand, with the increase of enclosure area in ecological construction projects, the number of livestock that has not been significantly reduced has increased grazing pressure outside the enclosure, resulting in some grassland degradation [46].

By comparing the evolution of NDVI and SIF, it was found that they had high consistency and spatial correlation. The simulation results agree that although the VPD in the source area has increased, the vegetation in most areas has become greener and more productive, consistent with previous research results [47]. However, due to the characteristics of the product, its change trend and the degree of influence by drought show certain regional differences. The results show that the impact of drought on SIF is slightly higher than that on NDVI, and the impact range is wider. In low vegetation coverage areas, NDVI is often significantly affected by soil background reflectance. Therefore, there is a certain lag and underestimation in monitoring the effects of drought on the activity of vegetation, and the SIF signal can monitor the impact of drought on vegetation activity more quickly and accurately than the vegetation index in terms of time and space changes. SIF can eliminate the influence of photosynthetically active radiation and is more sensitive to changes in water conditions. It can be used as an important indicator to monitor the occurrence and dynamic changes of regional drought [48,49].

4.3. Effects of Drought on Vegetation in the SRYY

Earlier research has indicated that SM is the most closely related terrestrial water variable to vegetation growth [50]. Consistent with the results of this study, low SM will seriously reduce vegetation greenness and weaken ecosystem carbon accumulation. In addition, VPD plays a crucial role in regulating the variation of vegetation in the source area. When most vegetation is exposed to high VPD stress, in order to reduce the risk of hydraulic failure triggered by intense transpiration, some or all of the stomata will be closed, eventually making the photosynthetic products insufficient [51]. It should be noted that high VPD stress will also increase the loss of water in moist soil, further aggravating the water depletion of plants and the risk of atmospheric drought.

Different vegetation types exhibit varying levels of adaptability to soil and atmosphere moisture. This study found that shrubs are the most sensitive to SM, and then sparse vegetation, grassland, and cultivated plants. Forest vegetation exhibits robust resistance to drought-induced stress. The response speed of grassland to SM was faster than that of forest land. This is due to the well-developed root system of the forest land, which can use the water stored in the deep soil under drought conditions to alleviate the effects of drought. Shallow root grasses usually absorb water from the upper soil layer, while shallow soil moisture is susceptible to drought, and the ability of the herbaceous plant xylem system to store water and carbon is low, resulting in poor resistance to drought stress [52]. In addition, in irrigated agricultural areas, efficient crop field management measures and irrigation techniques can improve crop water use efficiency [53]. In the meantime, vegetation provides positive feedback to soil and atmosphere moisture. The increase in vegetation can reduce soil evaporation and increase canopy area and density. Higher vegetation may increase soaking rain, stemflow and infiltration, thereby increasing root soil moisture [54].

In view of the dominant role of atmospheric drought and soil drought in vegetation growth, this study confirmed that SM in the SRYY had a stronger impact on vegetation greenness and productivity than VPD. Spatially, unlike other regions, the positive and negative causal effects of VPD on NDVI in the source area are equivalent, and the positive causal effect of forest land is significant. This could be attributed to the distinctive natural conditions within the alpine ecosystem The ecosystem in the SRYY is mainly limited by temperature, and the increase in VPD is primarily influenced by local temperature. As temperatures rise rapidly, the plateau forest land may have the capacity to adjust to the warming environment, lift the temperature limit in some areas and promote vegetation growth. At the same time, warming also causes more ice and snow to melt on the plateau, and the water available for plants to absorb is more sufficient, which indirectly makes the source vegetation less sensitive to the negative effects of atmospheric drought than other regions. With global warming, soil drought may still be the main water constraint limiting vegetation growth in the SRYY in the future, which is more apparent in the SRYA [20]. It is pointed out that in the past 40 years, the annual average temperature in the source SRYA the SRY has increased by about 0.8 °C and 0.7 °C, respectively, but the warming rate in the SRYA is still significantly greater than that in the SRY [55], which further confirms the results of this study.

4.4. The Response of Vegetation to Drought in Different Vegetation Types and Permafrost Degradation Areas

With the warming of the climate, the permafrost in the source area is gradually degraded. The phenomenon of permafrost degradation is usually associated with rising temperatures and freeze–thaw processes, which may lead to increased soil moisture. With the decrease of the frozen layer, the soil’s ability to store water is enhanced. Vegetation is sensitive to permafrost degradation [56]. This study pointed out that SM had a stronger effect on vegetation greenness and productivity than VPD, which was more obvious in the permafrost degradation area. On the one hand, in the frozen soil region, the freezing and thawing process of soil will form a frozen layer, which limits the storage of soil moisture and the acquisition of water by plants. However, permafrost degradation can lead to the melting of underground ice, and vegetation is relatively sensitive to this change when soil moisture increases [57]. On the other hand, vegetation in permafrost regions is typically constrained by water availability. The increase in soil moisture can provide a more suitable growth environment and promote the survival and reproduction of plants. Due to the limited soil moisture, vegetation may be more sensitive to the increase in water, which is reflected in the increase in growth rate, the adjustment of plant community structure and the increase in vegetation cover [58]. This response holds significant importance for the restoration of ecosystems in permafrost degradation areas and is of significance for the natural regeneration of vegetation, improvement of soil structure, mitigation of soil erosion, and improvement of regional ecological stability.

There is a strong interaction between frozen soil and vegetation in the source area. Studies have shown that in recent years, alpine meadow ecosystems are more sensitive to permafrost degradation than to climate change [59]. Since the late 1980s, the yellow-drying period of alpine meadows in the source area affected by permafrost degradation has been significantly delayed, resulting in a significant extension of the growth period and an overall decline in coverage. However, since 2004, the coverage increased and the aboveground biomass increased significantly [59], which was basically consistent with the results of this study. The water and heat state of the soil, that is, the temperature and humidity conditions of the soil, is an important link connecting the interaction between the two. On the one hand, in the source region and the surrounding permafrost regions, the observation data in the past 10 years show that the soil heat flux is on the rise, and the rate of change is from 1.0 W·m⁻²·(10a)⁻¹ to 2.9 W·m⁻²·(10a)⁻¹. With the increase in soil flux by 1.0 Wm⁻², the thickness of the permafrost active layer will increase by 24 cm [60]. The increase in soil heat flux may be the main direct reason for the decrease in freezing and thawing depth and freezing and thawing period of seasonal frozen soil in the plateau, and the increase in permafrost active layer and thawing period. On the other hand, due to the long freezing and thawing period of the frozen soil in the source area, with the increase of precipitation in the plateau in the past 60 years, precipitation mainly increases the temperature of the soil and accelerates the melting of the frozen soil [61]. Under the dual influence of global warming and human activities, the frozen soil and vegetation changes of SRYY are facing great challenges in the future, affecting the ecological security of the source area. Under various emission scenarios, it is shown that the future temperature in the source area will continue to rise and the precipitation will increase further [62]. Under this background, the soil temperature and humidity of the frozen soil in the source area show an increasing trend, and the permafrost area, the freezing depth and freezing period of the seasonal frozen soil will be further reduced, resulting in changes in the water cycle and water production in the basin. In addition, there is an obvious stratification in the soil at about 40~50 cm. The saturated zone in the upper layer is the main source of water absorption by plant roots. Under the background of increasing temperature, plant roots may develop to the deep soil layer, and vegetation NDVI will increase significantly [63,64]. This will have an important impact on the watershed vegetation and the local ecological environment.

5. Conclusions

This study reveals the temporal and spatial evolution of land–atmosphere water constraint events in the source region. The results show that the degree of atmospheric drought in the source region of the Yangtze River and the Yellow River increased from 2001 to 2021, and the soil drought was alleviated. The spatial and temporal distribution of NDVI and SIF is very consistent, indicating that SIF can better reflect the real state of vegetation growth. Both of them show the evolution process of ‘increase-decrease-increase’, and the vegetation coverage has improved significantly after 2008. Spatially, the vegetation coverage in the source area is characterized by high southeast and low northwest.

The PCMCI causal analysis method was also used to explore the effect of land–atmosphere water constraints on vegetation growth. The results showed that high VPD and low SM were water constraints limiting plant growth. Among them, shrubs are the most sensitive to soil drought, followed by sparse vegetation, grassland and cultivated plants. Forest vegetation has a strong ability to resist drought stress, and the response speed of grassland to SM was faster than that of forest land. The SM in the source area had a stronger effect on vegetation than VPD, and the permafrost degradation area was more obvious. SM dominated NDVI (SIF) changes in 62.87% (76.60%) of the permafrost degradation area.

Our research contributes to an enhanced understanding of the feedback relationship between global climate change and vegetation of different types. The study results provide scientific support for the scientific assessment and adaptation of permafrost, vegetation and climate change in the source area and provide reference for ecological protection in permafrost regions.