Spatio-Temporal Characteristics of Drought Events and Their Effects on Vegetation: A Case Study in Southern Tibet, China

Abstract

Frequent droughts in a warming climate tend to induce the degeneration of vegetation. Quantifying the response of vegetation to variations in drought events is therefore crucial for evaluating the potential impacts of climate change on ecosystems. In this study, the standardized precipitation index (SPI) was calculated using the precipitation data sourced from the China Meteorological Forcing Dataset (CMFD), and then the drought events in southern Tibet from 1982 to 2015 were identified based on the SPI index. The results showed that the frequency, severity, and intensity of drought events in southern Tibet decreased from 1982 to 2015, and the highest frequency of drought was found between 1993 and 2000. To evaluate the impact of drought events on vegetation, the vegetation characteristic indexes were developed based on the normalized difference vegetation index (NDVI) and the drought characteristics. The assessment of two drought events showed that the alpine grasslands and alpine meadows had high vegetation vulnerability (AI). The assessment of multiple drought events showed that responses of vegetation to drought were spatially heterogeneous, and the total explain rate of environmental factors to the variations in AI accounted for 40%. Among the many environmental factors investigated, the AI were higher at middle altitudes (2000–3000 m) than low altitudes (<2000 m) and high altitudes (3000–4500 m). Meanwhile, the silt soil fraction in the upper soil layer (0–30 cm) had the greatest positive correlation with AI, suggesting that areas with a high silt soil fraction were more sensitive to drought. The relative contribution rates of environmental factors were predicted by a multivariate linear regression (MLR) model. The silt soil fraction was found to make the greatest relative contribution (23.3%) to the changes in AI.

1. Introduction

Drought is considered to be one of the most complex and severe natural disasters affecting water resources, agricultural production, and natural ecosystems [1,2]. Water availability is the major driver of vegetation distribution and productivity in arid and semi-arid regions [3,4]. Previous studies have suggested that arid/semi-arid ecosystems dominate the inter-annual variation in terrestrial CO₂ sinks and are highly sensitive to drought trends [5,6,7]. Therefore, understanding the responses of terrestrial ecosystems to drought is of great importance for assessing the vulnerability of vegetation to extreme climate events and for ensuring the stability of terrestrial ecosystems.

The identification and characterization of drought events is crucial to analyze the responses of vegetation to drought [8]. Previously, several methods have been used to characterize drought including the threshold method [9], the runs theory [10,11], the empirical orthogonal method [12,13], the wavelet analysis [14], and principal component analysis [15]. These methods usually quantify the characteristics of drought based on the time evolution in a fixed region or the spatial pattern over a fixed duration. This could limit the comprehensive understanding the spatio-temporal structure of drought events [16]. To overcome this limitation, a three-dimensional clustering method was developed [17]. It defines a drought event as a continuous cluster in the three-dimensional space in which the drought index is below a specified threshold. The application of this method could improve our understanding of the development of drought events, the changes in drought characteristics, and the potential effects of drought on vegetation health.

The selection of an appropriate drought index is the basis of the identification and characterization of drought [18]. The most frequently used indices include the standard precipitation index (SPI) [19], the standard precipitation evapotranspiration index (SPEI) [20], and the Palmer drought severity index (PDSI) [21]. The SPI is calculated based on precipitation, and it reflects drought conditions over different timescales. In addition to precipitation, SPEI also considers the effects of temperature on drought condition. Uncertainties in the SPEI index were controlled by the potential evapotranspiration data [22]. The PDSI was originally developed to measure cumulative deviations from the surface water balance. However, the PDSI cannot identify drought over different timescales because it lacks the multi-scalar characteristic of droughts. Compared with SPEI and PDSI, the SPI is simpler, multi-scalar, and less affected by the uncertainty of input data [22].

Among the various vegetation indices, the normalized difference vegetation index (NDVI) is an efficient indicator for monitoring vegetation conditions (e.g., land degradation and desertification) [23] and is frequently used in drought detection and assessment [24,25]. For example, Vicente-Serrano [26] used the NDVI to study the response of global land biomes to drought timescales. Ji and Peters [27] studied the response of vegetation to drought conditions using NDVI data in the American Great Plains. Over the past decades, the NDVI index has been widely used to assess the relationships between vegetation and drought conditions. Recently, drought has become a temporary feature of the climate. Michele [28] noted the NDVI anomaly could be used to detect anomalous vegetation changes, compared with the assessment based on the average or “normal” situation. Therefore, the NDVI anomaly is frequently used to assess the impacts of climate events on the vegetation situation.

The ecosystem of southern Tibet is fragile due to its specific geological, topographic, and climate conditions [29]. Thus, southern Tibet is often influenced by droughts, storms, frost, and other natural disasters. Frequent droughts would further aggravate the shortage of water resources and have negative impacts on vegetation growth [30]. Observations from southern Tibetan meteorological stations have indicated that the occurrence of drought events has generally increased in the humid southern part of South Tibet and decreased in the arid part of southern Tibet [31]. However, in southern Tibet, our understanding of the responses of vegetation to drought remain unclear on a large scale [32]. Studying the response of vegetation to drought events will improve our understanding of the land–air interaction processes in the earth system model which would be beneficial for ecosystem stability.

The objective of the study was to investigate the impacts of drought on vegetation through the analysis of the relationship between drought characteristics and vegetation index. The main works were focused on three major areas: (i) the temporal and spatial characteristics of drought events in southern Tibet; (ii) the changes in responses of vegetation to drought severity for different vegetation types and elevation bands through comparing two drought events; (iii) the impacts of exogenous factors on vegetation vulnerability. The results furthered our understanding of temporal and spatial variation of drought characteristics and the potential impacts of drought events on vegetation.

2. Materials and Methods

2.1. Data Sources

2.1.1. Study Area

The study area (26°49′N–33°47′N and 78°23′E–97°34′E) was located in the southern part of the Tibetan Plateau with a length of 1700 km from east to west and a width of 1000 km from south to north. This area is bounded by the Indian Himachal Pradesh in the west, the Hengduan Mountains in the east, the Gangdise and Nyenchentanglha Mountains in the north, and the Himalayas in the south (Figure 1). The average elevation of the study area is more than 4000 m, and the total area is 520,000 km².

The southeastern part of the study area is the hottest area of the Tibetan Plateau with an annual average temperature of 10–25 °C and an annual rainfall of 1000–3000 mm. The central part (Lhasa to Shigatse) is within a temperate and semi-arid climate zone in which more than half of the region has an annual average rainfall of 400–600 mm and an annual average temperature of −5 to 10 °C. The northwest part is in the plateau sub-frigid climate zone with annual rainfall of less than 400 mm and an annual average temperature of −5 to 0 °C [33]. The interaction between water and heat plays a decisive role in the vegetation distribution. From the southeast to the northwest, the vegetation types change from subtropical coniferous forests and broad-leaved forests to alpine grassland and desert vegetation.

2.1.2. Meteorological Data

The China Meteorological Forcing Dataset (CMFD) was the first near-surface meteorological dataset developed for land surface processes in China with the advantage of high spatial–temporal resolution. This dataset integrates five auxiliary data sources, including in situ station data, tropical rainfall observation mission (TRMM) precipitation data, global energy and water-cycle experiment project (GEWEX-SRB) downward short-wave radiation data, Princeton meteorological forcing data, and global land data assimilation system (GLDAS) data. Its record spans the period from 1979 to 2018 with a temporal resolution of 3 h and a spatial resolution of 9 × 9 km (

Table 1. Environmental variables used in the study.

Classification and Data Source	Temporal Resolution	Spatial Resolution	Description
Climate variables westdc.westgis.ac.cn/data	3 h	9 km	Precipitation
Vegetation indexes http://ecocast.arc.nasa.gov/data/pub/	15 d	8 km	NDVI (normalized difference vegetation index)
Soil variables https://data.tpdc.ac.cn/zh-hans/	-	1 km	Soil texture (i.e., sandy soil, silt soil, clay soil and gravel)
Topographic variables http://www.gscloud.cn	-	90 m	Altitude, relief, and aspect
Vegetation types http://www.resdc.cn	- (vector)	- (vector)	Basic vegetation types

). The CMFD contains all seven near-surface meteorological elements applied for land modeling, including temperature, air pressure, air specific humidity, total wind speed, downward short-wave radiation, downward long-wave radiation, and precipitation rate. This dataset has been widely used in plant productivity estimation studies [34], vegetation growth analyses [35], and lake area simulations [36]. In this study, the 3 h precipitation data were processed into monthly total precipitation.

2.1.3. Global Inventory Modelling and Mapping Studies (GIMMS) NDVI

The NDVI is an indicator of canopy greenness and has been used to monitor vegetation dynamics at regional and global scales. Calculated from spectral reflectance in visible/near infrared wavebands, the NDVI has proved to be sensitive to changes in biomass, canopy structure, and water content, etc. [37,38,39]. This GIMMS NDVI data has an 8 km resolution and a bi-weekly temporal resolution covering the period of 1981–2015 (Table 1). The data has a longer time series and a wider spatial range than the Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI and Satellite Pour l’Observation de la Terre’ VEGETATION (SPOT/VGT) NDVI. Artifacts caused by orbital and solar and zenith angle changes have been calibrated by the empirical mode decomposition/reconstruction method [40]. Uncertainties caused by volcanic aerosol has been eliminated from the GIMMS NDVI using atmospheric corrections [41]. The GIMMS NDVI data were used to generate monthly NDVI data using the maximum value composite method, and the data were resampled to 9 km using the nearest neighborhood method. Pixels with an annual average NDVI of less than 0.1 were considered to be non-vegetated pixels and were not included in the subsequent calculations [42].

2.1.4. Topographic Parameters

The Shuttle Radar Topography Mission digital elevation model (SRTM DEM) was constructed using radar interferometry (InSAR) over 10 days in February 2000. The SRTM DEM was provided freely by the China Geo-Spatial Data Cloud (http://www.gscloud.cn) with a resolution of 90 m (Table 1). Three digital terrain variables were obtained by SAGA GIS (version 2.3.0, System for Automated Geo-Scientific Analyses, http://www.saga-gis.org) based on the SRTM DEM including relief and aspect.

The soil moisture content of the shady slopes was higher than that of sunny slopes, and the soil moisture of half shady and half sunny slopes had values between the two extremes. We classified the slope into four types (i.e., sunny slope, shady slope, half sunny slope, half shady slope) [43,44] in which 135–225° was a sunny slope, 225–315° was a half sunny slope, and 0–45° and 315–360° were the shady slopes, while the remaining slopes were half shady slopes. Topographic relief was mapped by measuring the surface slope within a 3 × 3 pixel moving window. A fishnet with a spatial resolution of 9 × 9 km extracted the average altitude, average relief, and aspect based on the mean function of zonal statistics from the ArcGIS software. The fishnet was converted to a raster with a spatial resolution of 9 km.

2.1.5. Soil Texture

Data on soil texture were sourced from the soil properties dataset for land surface modeling over China provided by the National Tibetan Plateau/Third Pole Environment Data Center (https://data.tpdc.ac.cn/zh-hans/) (Table 1). Soil texture affects the distribution and duration of water stored in soil, both of which affect plant distribution [45,46]. The gridded soil texture had a spatial resolution of 1 km. Soil information was obtained for eight standard layers. The soil profiles were divided into eight standard layers (<4.5 cm, 4.5–9.1 cm, 9.1–16.6 cm, 16.6–28.9 cm, 28.9–49.3 cm, 49.3–82.9 cm, 82.9–138.3 cm, and 138.3–229.6 cm). The data for the top layers (0–0.289 m) were used. The fishnet with a spatial resolution of 9 × 9 km was used to extract the average percentages of sand soil, silt soil, and clay soil based on the mean function of zonal statistics from the ArcGIS software. Soil texture (percentage of sand soil, silt soil, and clay soil), as the static variable, was used to explore the effects of soil factors on vegetation vulnerability under drought conditions.

2.1.6. Vegetation Type

The vegetation types used in this study were taken from the 1:1,000,000 vegetation type data provided by the Resource and Environmental Science and Data Center (http://www.resdc.cn) (Table 1). The vegetation was divided into 11 different types, namely, grassland, meadow, desert, shrub, alpine vegetation, broad-leaved forest, coniferous forest, cultivated vegetation, alpine vegetation, and coniferous and broad-leaved mixed forest.

2.2. Methodology

2.2.1. Calculation of the SPI

The impact of drought is difficult to accurately identify and quantify because of the complexity of the process. Therefore, some drought indicators, such as rainfall, soil moisture, and runoff, are often used to assess drought conditions. The commonly used indicators include the SPI, SPEI, and standardized soil moisture Index (SSMI) [47,48,49]. The SPI is a reference drought index recommended by the World Meteorological Organization [50]. It is determined using the standardized precipitation anomaly under a given probability distribution and can be calculated on multiple time scales (usually 1, 3, 6, 12 or 24 months) [51]. It has been shown that SPI-3 has the greatest correlation with vegetation change and drought evolution. The formula for the SPI index is as follows:

$$\mathrm{f}\left({\mathrm{x}}_{\mathrm{k}}\right)=\frac{1}{{\mathsf{\beta }}^{\mathsf{\alpha }}\Gamma \left(\mathsf{\alpha }\right)}{\mathrm{x}}_{\mathrm{k}}^{\mathsf{\alpha }-1}{\mathrm{e}}^{-\mathrm{x}/\mathsf{\beta }}$$

(1)

where $\mathsf{\alpha }$ is the form parameter ($\mathsf{\alpha }$ > 0), $\mathsf{\beta }$ is the scale parameter ($\mathsf{\beta }$ > 0), ${\mathrm{X}}_{\mathrm{k}}$ is the cumulative rainfall for consecutive months (mm), $\Gamma \left(\mathsf{\alpha }\right)$ is the Gamma function, and $\overline{\mathrm{x}}$ is the annual average rainfall amount; $\mathsf{\beta }$ and $\mathsf{\alpha }$ are calculated by maximum likelihood methods.

$$\mathsf{\alpha }=\frac{1+\sqrt{1+4\mathrm{A}/3}}{4\mathrm{A}}$$

(2)

$$\overline{\mathsf{\beta }}=\overline{\mathrm{x}}/\mathsf{\alpha }$$

(3)

$$\mathrm{A}=\ln \left({\overline{\mathsf{x}}}_{\mathrm{k}}\right)-\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}\ln \left({\mathrm{x}}_{\mathrm{k}}\right)$$

(4)

According to Equation (5), the cumulative probability under a given time scale is calculated as follows:

$$\mathrm{F}\left(\mathrm{x}\right)={\int }_0^{\mathrm{x}}\mathrm{f}\left({\mathrm{x}}_{\mathrm{k}}\right)\mathrm{dx}=\frac{1}{\Gamma \left(\mathsf{\alpha }\right){\mathsf{\beta }}^{\mathsf{\alpha }}}{\int }_0^{\mathrm{x}}{\mathrm{x}}^{\mathsf{\alpha }-1}{\mathrm{e}}^{-\mathrm{x}/\mathsf{\beta }}\mathrm{dx}$$

(5)

The “SCI” package in the R statistical software [52] was used for the calculation of SPI. The classification of SPI is shown in

Table 2. Classification of the standardized precipitation index (SPI).

SPI Intervals	SPI Classes
SPI > 2	Extremely wet
2 > SPI > 1.5	Severely wet
1.5 > SPI > 1	Moderately wet
1 > SPI > −1	Normal
−1 > SPI > −1.5	Moderately dry
−1.5 > SPI > −2	Severity dry
SPI < −2	Extremely dry

. The SPI value of −1 has been used as the threshold for extracting drought events in many regions, such as worldwide [53], the United States [17], Europe [16], Central Asia [54], China [55], the Pearl River Basin [56], Southwest China [57], and Poyang Lake Basin [58].

2.2.2. Drought Identification

A three-dimensional clustering algorithm was used to identify the drought object in three-dimensional space (longitude, latitude, and time). The three-dimensional space had a size of V_lat × V_lon × V_t, where V_lat and V_lon represented the grid numbers along with the latitude and longitude, while V_t stood for the month number during the period of 1982–2015.

The three-dimensional cluster algorithm is frequently used to identify the non-overlapping and independent extreme climate events on time and space scales. In this study, the drought index maps were first converted to binary distribution maps according to the threshold of drought intensity, in which pixels with an SPI index less than the threshold of −1 threshold were assigned the value 1, and pixels greater than −1 were assigned the value 0. After that, the connected domain of the monthly binary distribution map was analyzed, and the adjacent drought patches with a value of 1 were identified. The patches smaller than the fixed threshold were filtered out because they may affect the identification of large-scale drought events. The threshold area of approximately 8000 km² was determined as 1.6% of the total study area of the study region, which was about 8000 km² according to previous studies [8,59]. The overlapping area between the drought patches for adjacent months was calculated and a threshold area of 6400 km² for the overlapping area was specified [53]. When the overlapping area was larger than 6400 km², the drought patches for adjacent months were marked as the same drought event. Otherwise, they were considered to be different drought events. Following the above procedure, all the drought patches appearing in adjacent months were traversed until all the drought patches were marked, and different marks represented different drought events.

2.2.3. Drought Characteristics

The drought characteristics can be calculated using a descriptive statistical algorithm such as the occurrence time, end time, duration (DD), area (DA), centroid (DC), severity (DS), and intensity (DI). The first five indexes were generated by the regionprops function in MATLAB, but the drought severity (DS) and intensity (DI) indexes needed to be further calculated.

(1) Drought severity refers to the cumulative value of the drought index from the first month to the last month of the drought event and mainly reflects the degree of water scarcity. The drought severity of the nth drought event can be calculated according to Equations (6) and (7).

$$D{S}_n={\displaystyle \sum _i^{N_{lat}}\sum _j^{N_{lon}}\sum _k^{N_t}}s\left(i,j,t\right)$$

(6)

$$\begin{array}{l}s\left(i,j,k\right)\\ =\{\begin{array}{ll}A{I}_{i,j,k}\times ARE{A}_n\times 1month,& ifgrid\left(i,j\right)islabelledbythenthdroughtevent\\ 0,& ifgrid\left(i,j\right)isnotlabelledbythenthdroughtevent\end{array}\end{array}$$

(7)

where i and j represent the horizontal and vertical coordinates of the corresponding pixel, K represents the time coordinate of the corresponding pixel, and AI_i,j,k represents the SPI drought index at the corresponding pixel position.

(2) Drought intensity refers to the average drought severity per unit time and the area of the drought event and can be calculated from Equation (8).

$$DI=D{S}_n/\left(D{A}_n\times D{D}_n\right)$$

(8)

where DD_n represents the number of months of the nth drought event, and DA_n represents the areas affected by the nth drought event (km²).

2.2.4. Vegetation Anomaly Index

The vegetation anomaly is based on the NDVI time series and can reflect vegetation health over a short timescale. The trends and seasonality of an NDVI time series can affect the calculation of the vegetation anomaly value. First, the empirical mode decomposition (EMD) was used to remove the inter-annual trend components with a period of more than 3 years in the NDVI time series. Then, the Z-score standardized method was used to calculate the monthly NDVI anomaly from 1982 to 2015, because it effectively diminishes the influences of seasonal changes.

$$NDV{I}_{anomaly,i,j}=\frac{NDVI\left(i.j\right)-\overline{NDV{I}_i}}{{\delta }_i}$$

(9)

where $NDV{I}_{anomaly,i,j}$ is the standard deviation of the NDVI value for the jth month during the ith year, $\overline{NDV{I}_i}$ is the averaged NDVI for the ith month over all years, and ${\delta }_i$ is the variance of the NDVI sequence for the ith month over a number of years. Anomalies were classified into positive anomalies (NDVI_anomaly > 1), normal (|NDVI_anomaly| < 1), and negative anomalies (NDVI_anomaly < −1). Negative NDVI anomalies indicated that vegetation growth was negatively affected.

2.2.5. Vegetation Characteristics under Drought Stress

Three vegetation characteristics indexes were established under drought stress in order to effectively evaluate the negative effects of drought events on vegetation changes. They are the area percent of the negative vegetation anomaly (AP), the negative vegetation anomaly severity (AS), and vegetation vulnerability (AI). The first two indicators were used to evaluate the characteristics of the vegetation negative anomaly under a single drought event, and the third indicator was used to evaluate the long-term possibility of a negative vegetation anomaly. The specific descriptions are as following:

(1) Area percent of the negative vegetation anomaly (AP) indicates the area of the negative vegetation anomaly out of the total areas affected by drought. The AP index was calculated using Equation (10).

$$A{P}_i=\frac{N{T}_i}{T{P}_i}$$

(10)

where NT_i represents the number of negative anomaly pixels in the ith drought event, and TP_i represents the total number of pixels affected by the ith drought event.

(2) The negative vegetation anomaly severity (AS) indicates the degree to which the vegetation is affected by a water storage deficit. The AS index was calculated using Equations (11) and (12).

$$A{S}_n={\displaystyle \sum _i^{N_{lat}}\sum _j^{N_{lon}}\sum _k^{N_t}}NS\left(i,j,k\right)$$

(11)

$$\begin{array}{c}NS\hfill \\ =\{\begin{array}{ll}NDV{I}_{anomaly,i,j}\times area\times 1month,& ifgrid\left(i,j\right)islabelledbynthdroughtevent\\ 0,& ifgrid\left(i,j\right)isnotlabelledbynthdroughtevent\end{array}\hfill \end{array}$$

(12)

where n is the label of a drought event, area represents the grid area, 1 month is the time period, and NS_n is the monthly negative vegetation anomaly severity of the nth drought event.

(3) Vegetation vulnerability (AI) indicates the possibility of negative anomaly of vegetation at a specific location for a certain period under drought stress. The AI index can be calculated as Equation (13).

$$A{I}_{i,j}=\frac{N{F}_{i,j}}{T{F}_{i,j}}$$

(13)

where TF_i,j indicates the total number of drought events at pixel (i,j), and NF_i,j indicates the number of negative vegetation anomalies at pixel (i,j).

2.2.6. Multivariate Linear Regression (MLR)

Multivariate linear regression is widely used for analyzing the correlations between the variation in ecosystem and the variation in environmental factors [60,61]. In a multivariate regression analysis, the assumptions of the regression analysis are normal distribution, linearity, and no multicollinearity among the independent variables. Multivariate regression analysis was formulated as follows:

$$\mathrm{AI}={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_n{x}_n+\epsilon$$

(14)

where AI represents the vegetation vulnerability under drought stress, $x_n$ represents the environmental variables, ${\beta }_n$ represents standard regression coefficients, and $\epsilon$ represents the constant of the regression model. In the multiple regression model, the p-value for an F-test determines the addition and removal of terms [62]. For an MLR, all the variables were chosen to model the correlation between AI and environmental factors. According to the standard regression coefficients of the MLR, the relative contribution rate (PR_i) of each variable was expressed as:

$$P{R}_i=\frac{\left|{\beta }_i\right|}{{\sum }_{i=1}^n\left|{\beta }_i\right|}$$

(15)

2.2.7. The Overall Flowchart of the Methodology

The framework of the drought identification and analysis methodology is shown in Figure 2. The first process, including step 1 to step 3, was used to calculate the drought characteristic indexes. The second process (i.e., step 4) was to estimate the vegetation characteristics under drought stress.

3. Result

3.1. Temporal Characteristics of Drought Events

A total of 44 drought events were identified in southern Tibet between 1982 and 2015. These drought events lasted 3, 4, 5, 6, 7, 9, and 11 months, respectively, with the corresponding frequencies of 18, 5, 7, 5, 4, 2, and 2 times. Most of the drought events lasted from 3 to 5 months, and only two lasted for 11 months. The longest two drought events occurred between January and November in 1998 and from January to November in 2004.

The changes in drought characteristics over the lifetime of a drought (i.e., DD, DS, and DI) and for each monthly step (i.e., AP, AS, and AI) are plotted in Figure 3. The DS, DA, and DI are represented by solid black rectangles, while MS, MA, and MI are represented by the colored cylinders. From 1982 to 2015, the severity and scope of the drought events showed a downward trend, with tendency rates of 1.85 × 10⁴ year⁻¹ and 4192.42 km²·year⁻¹, respectively. Furthermore, the drought intensity fluctuated around the average state without a consistent trend. The drought frequency from 1982 to 2015 was 1.26 times·year⁻¹, with the highest frequency being 1.86 times·year⁻¹ between 1993 and 2000. After 2000, short-term and small-scale drought events were dominant. Short-term (3 months) drought events often started in spring and summer and at the end of the current season. Mid-term drought events (4–7 months) mostly started in summer and ended in autumn and winter. Long-term drought events (9–11 months) started in autumn and winter (Figure 4).

3.2. Spatial Characteristics of Drought Events

The spatial distribution of the centroids of the 44 drought events are shown in Figure 5a. The size of the circle represents the drought severity, and the color represents the duration of the drought. Figure 5b depicts the spatial distribution of the drought events. The beginning of the arrow indicates the location of the weighted centroid based on severity during the first month of the drought event, and the tail of the arrow indicates the location of the weighted centroid based on the severity during the last month of the drought event. The development of drought events was dominated by the horizontal direction, which is independent of the duration and probably related to the channels for water vapor transport. There are many meridional canyons over the southern Tibetan Plateau (STP) that can serve as channels for water vapor transport [63].

3.3. Evaluation of a Single Drought Event

3.3.1. Effects of a Single Drought Event on Vegetation

The characteristics of drought events were determined in the previous section. The amplitude and spatial extent of the drought event in 1982 exceeded any other event, affecting more than 70% of the area of southern Tibet. The core of the drought event in 1982 was mainly in the central part of southern Tibet, whereas the western and eastern parts of southern Tibet were less impacted, exhibiting a lower drought severity (Figure 6a). Figure 6b displays the spatial pattern of the drought severity in 1987. The areas affected by the drought event in 1987 was similar to those affected in 1982, but the average severity (252 km²·year) was only half of that in 1982. Figure 6c showed the difference in severity of the two drought events. The positive values of the R-DS suggest that the severity of the drought event in 1982 exceeded that in 1987 and vice versa. As shown in Figure 6c, the areas with positive values covered most of the area of Shigatse, Shannan, and Lhasa, accounting for over 70% of the total studied area, while the areas with negative values were only concentrated in the Ali region in the west.

The R-NDVI calculated from the subtraction method using the average growing season (May–October) NDVI in 1982 and 1987 represent the differences in vegetation greenness of the two events (Figure 7a,b). As displayed in Figure 7c and Figure 8b, the positive areas were located in the central and western parts, accounting for 55.7% of the total studied area, suggesting that the average growing-season NDVI of these areas in 1982 was lower than that in 1987.

3.3.2. Different Responses of Vegetation Types and Elevation Bands

The response of vegetation types to water deficit is a crucial issue for the geographical pattern of vegetation and is a central concept to understanding the structure and dynamics of terrestrial ecosystems [26,64]. Based on the relationships between R-DS and R-NDVI of different vegetation types, we found that the mean and median of the R-NDVI of alpine grassland increased with the increase in the R-DS, suggesting that alpine grassland was more sensitive to drought (Figure 9a). As the R-DS value increased, the range of the R-NDVI of alpine meadows became stable, maintaining a value between 0–0.02 (Figure 9b). Unlike alpine grassland and alpine meadow, there was a threshold for the R-DS in sub-alpine broad-leaved shrubs (Figure 9c). When the R-DS was higher than 800 km²·month, most of the R-NDVI values were above zero, suggesting that vegetation greenness significantly decreased. The R-NDVI for broad-leaved forest and coniferous forest fluctuated around zero (Figure 9e–f), which was not consistent with the R-DS values, suggesting that there was little influence of drought on vegetation activity. This could be related to the physiological characteristics of woody plants, which allows them to reduce the damage caused by water deficits [65].

The relationships between R-NDVI and R-DS for different elevation bands are shown in Figure 9. When the elevation is lower than 4000 m, the R-NDVI values fluctuated around zero as the R-DS increased, suggesting that the negative impacts of drought might be limited (Figure 10a–c). When the elevation is between 4000–5000 m, together with the R-DS exceeding 400 km²·month, vegetation greenness decreased significantly (Figure 10d,e). However, when the elevation was higher than 5000 m, the R-NDVI data for most vegetation pixels were higher than 0, suggesting that vegetation in high elevation areas was more sensitive to drought (Figure 10f), which was consistent with some other studies [66,67].

3.4. Assessment of the Effects of Multiple Drought Events Effects on Vegetation

3.4.1. Relationships between Drought and Vegetation Characteristics

The relationships between the drought characteristics (DD, DA, and DS) and the vegetation characteristics (AP and AS) were analyzed in order to better understand the response patterns of vegetation to drought. Figure 11 shows a weak correlation between AP and drought characteristics (R² < 0.1, p-value > 0.01). This suggested that vegetation in some areas was relatively resistant to drought stress. However, the drought characteristics can explain a large amount of the AS (between 20% and 50%, p-value < 0.01), suggesting that vegetation in some areas sensitive to drought would decline as the drought severity increased.

The above differences were dependent on two factors: (1) the different terrain, soil, and long-term climatic conditions result in the different responses of vegetation to drought, and (2) the water demands of the different vegetation types were different. In summary, environmental conditions may also promote or mitigate the negative effects of water deficits on vegetation.

3.4.2. Effects of Environmental Factors on the Vegetation Response to Drought Stress

An MLR model was used to determine the effects of environmental factors (such as soil texture, vegetation types, and topography) to the spatial variation in vegetation vulnerability (AI). The results revealed that all independent variables were significant statistically (p < 0.01). Meantime, all of the independent variables accounted for 40% of total variation in AI as a dependent variable, with no multi-collinearity due to the variance inflation factor (VIF) values varying from 1 to 5 (

Table 3. Estimated parameter and significance levels in the multivariate linear regression analysis. Variance inflation factor (VIF).

Model	Coefficients		t-Test	p-Value	VIF
	B	SE
Constant	37	4.098	9.249	3.36 × 10−20
Elevation	−0.005	0.000	−11.808	9.85 × 10−32	5.772
Relief	−0.142	0.022	−6.464	1.12 × 10−10	2.732
Sunny slope	0.198	0.042	4.694	2.75 × 10−6	1.354
Shady slope	−0.175	0.053	−3.28	1.00 × 10−3	1.4
Silt	0.68	0.08	8.397	6.00 × 10−17	2.035
Gravel	−0.216	0.045	−4.845	1.30 × 10−6	1.768
Alpine grassland	−0.118	0.01	−12.27	4.15 × 10−34	1.94
Temperate desert grassland	−0.28	0.016	−18.05	1.73 × 10−70	1.239
Temperate desert shrub	−0.316	0.024	−13.4	3.00 × 10−40	1.17
Alpine sparse vegetation	−0.139	0.011	−12.65	4.03 × 10−36	2.45
Deciduous broad-leaved shrub	−0.128	0.023	−5.55	3.00 × 10−8	1.085
Cultivated vegetation	−0.202	0.038	−5.327	1.04 × 10−7	1.41
Evergreen broad-leaved shrub	−0.044	0.011	−3.95	8.00 × 10−5	1.457
Temperate typical grassland	0.1	0.024	4.14	3.49 × 10−5	1.19
Broad-leaved forest	−0.176	0.027	−6.55	6.23 × 10−11	2.206

).

To further understand the impact of environmental factors on AI, the relationships between different environmental factors and AI were also investigated at the pixel scale. Figure 12a showed that change in altitude had a large impact on AI. When the altitude was lower than 2000 m, it had a positive correlation with AI, while the correlation was significantly negative when the elevation was between 3000 m and 4500 m (Figure 12a). This result suggested that the AI of vegetation in middle elevation areas (2000–3000 m) might be higher than that of low elevation areas (<2000 m) and high elevation areas (3000–4500 m). When the topographic relief was less than 1500 m, there was a positive correlation between topographic relief and AI, but when the topographic relief increased continually, the AI may have decreased (Figure 12b). In terms of slope aspect, the correlations between the area fractions of both the sunny and shady slopes and AI were insignificant (Figure 12c,d).

Relationships between soil texture and AI were shown in Figure 13 and negative correlation between the volume fraction of gravel and the AI was found. The larger the content of coarse particles, the smaller the AI, indicating they were less sensitive to drought. In contrast, the higher the volume fraction of silt, the greater the AI, suggesting a stronger sensitivity to drought. Walter et al. [68] found that sandy soils are beneficial to the growth of vegetation in arid regions. Noy-Meir [69] used the inverse texture effect to explain the phenomenon of “more dense perennial vegetation” in coarse soil.

3.4.3. Relative Contributions of Environmental Factors to AI

The relative contributions made by fifteen important environmental factors to the AI were assessed using the standard regression coefficients from the MLR (Figure 14). The results showed that silt content, sunny slope, and temperate typical grassland had positive effects on the AI, and their relative contribution rates were 23.3%, 6.8%, and 3.4%, respectively (Figure 15). The rest factors had negative effects on the AI, and the area fractions of temperate desert grassland and temperate desert shrubs accounted for 9.6% and 11% of the variations in AI, respectively. The relatively contribution of silt content to AI was the highest. This may be due to the high silt content of the surface soil, leading to an increase in water evaporation and a rise in drought severity.

4. Discussion

4.1. Regional Climate Characteristics

The results outlined in Section 3.1 showed that the frequency, range, and severity of drought events had a significant downward trend at the inter-annual scale. This suggested that southern Tibet is undergoing a wetting process, which is inconsistent with the overall warming and drying trends in China. Some studies have suggested that regional climate variability may be influenced by atmospheric circulation patterns and the elongation storm track. Zhu et al. [70] suggested that the extreme dryness in the Tibetan Plateau was associated with an anomalous cyclone over South and Southeast Asia, and the extreme wetness was affected by a zonal elongation of the North Atlantic storm track. Wang et al. [71] found that the Tibetan Plateau is experiencing a significant wet trend, at both the annual and seasonal scale, that is affected by enhancements to the summer Indian Ocean monsoon via the analysis of the summer drought time series. However, contradictory results were also reported based on the precipitation data recorded in the Indian Ocean monsoon region during the most recent 60 years, probably due to the weakening intensity of the Indian Ocean monsoon [72]. These different results may be caused by complex geomorphological features. Liu and Yin [73] analyzed the relationship between inter-annual precipitation and atmospheric oscillations in different regions of the Tibetan Plateau and concluded that the North Atlantic Oscillation may be the major reason for the precipitation differences in precipitation between the northern and the southern parts of the Qinghai Tibet Plateau at the inter-decadal scale. Since the mid-nineteenth century, the North Atlantic Oscillation and the Pacific inter-decadal Oscillation have been in a negative phase, which has amplified the impact of the La Niña phenomenon (El Niño Southern Oscillation (ENSO) negative phase) on the Indian Ocean monsoon. This increased the transport capacity of water vapor from the Bay of Bengal, resulting in a decreasing trend in the frequency of drought events in southern Tibet.

4.2. Effects of Environmental Factors on Vegetation Responses to Drought Stress

The local environment also plays an important role in vegetation evolution, and often leads to ecotype differentiation within species and changes the responses of vegetation to drought. The altitude gradients were considered to be one of important factors affecting the vegetation vulnerability [74]. Some studies suggested that the altitude gradients affect the vegetation vulnerability by altering the vegetation cover, and the larger the vegetation cover, the smaller the vegetation vulnerability [75,76]. Therefore, when the altitude increased, vegetation cover decreased, and the vegetation vulnerability strengthened. In addition, some other studies suggested that the positive relationship between vegetation vulnerability and altitude gradient might be attributable to the water-use efficiency [77,78]. The water availability and air temperature declined with the altitude increase, leading to a relatively slow recovery rate of vegetation after drought. However, most studies were carried out in the region with the average altitude below 3000 m. The result outlined in Section 3.4.2. showed that the vegetation vulnerability might be higher in the middle altitude region (2000–3000 m) than low altitude region (<2000 m) and high-altitude region (3000–4500 m). This result might be related with vegetation types and air temperature. The basic vegetation types in low-altitude areas (<2000 m) over the southern Tibet are broad-leaved forest and coniferous forest, which are more resistant to drought due to the fact of their physiological structure. When the altitude exceeded a certain height, the impact of thermal factors on vegetation change was prominent. Lower temperature leads to a slower metabolism rates, resulting in the accumulation of large amounts of non-structural carbohydrates in the vegetation. These carbohydrates (including fructose, glucose, and fats) can maintain plant growth and regulate osmotic pressure in cells under drought stress [79,80].

The other two environmental factors closely related to the changes in drought sensitivity were soil texture and vegetation types. Soils with high silt content can store more water near the surface than soils with coarser texture, and then be subjected to high evaporative losses, resulting in a relatively transient water resource in the topmost soil layers. Conversely, the high percentage of porosity in the coarse topsoil allows the rainwater to percolate rapidly downwards. Water in deep soil provides more opportunities for vegetation growth [81]. Duan et al. [82] have used a hydrogen and oxygen isotope tracer method to study the changes in water resources in the topmost soil and subsoil layers of forest, shrub, and desert ecosystems across the Tibetan Plateau. The results showed that during the dry season, the ¹⁸O isotope ratio in the topmost soil layers changes considerably, because of the strong evaporation. While in subsoil layers, the variation of the ¹⁸O isotope ratio was relatively small due to the stability of water storage in subsoil layers. Therefore, together with the infiltration of rainwater and the water stability in deep soil, soil with high contents of coarse-grained particles usually leads to a high-water storage level in the deep soil. This can increase the buffering capacity of local vegetation to drought.

Different vegetation types have different physiological mechanisms, which are key factors enabling different vegetation to adapt to changing water availability. The results outlined in Section 3.3.2. showed that alpine grassland and alpine meadow were more sensitive to drought severity. These results were consistent with previous studies. They reported that grassland and meadow areas were more sensitive to drought and had functional strategies enabling them to adapt to the changing water availability [83,84,85]. Baldocchi et al. [86] and Wolf et al. [87] reported that forests had a higher drought resistance than grasslands, due to the endogenous factors (i.e., deep root systems). Meanwhile, drought sensitivity among vegetation types is also driven by exogenous factors (such as topography and soil) [88].

5. Conclusions

This study evaluated the space–time characteristics of drought events in southern Tibet based on a three-dimensional clustering algorithm. The SPI drought index was used to construct the three-dimensional space of longitude × latitude × time. Several drought characteristics, such as duration (DD), area (DA), severity (DS), and intensity (DI) were used to depict the spatio-temporal variation. To evaluate the effects of drought on vegetation, three vegetation characteristic indicators (such as AP, AS, and AI) were used. The main conclusions are as follows:

(1) At the inter-annual time scale, the frequency of drought events showed a decreasing trend, which may be related with the La Niña phenomenon since the mid-nineteenth century. Seasonally, more drought months occurred in summer and autumn. Spatially, the direction in which the drought events developed over the study area was dominated by an east–west trajectory because of the channels for water vapor transport;

(2) Based on the comparison of two drought events in 1982 and 1987, the alpine grasslands and alpine meadows were more vulnerable to drought severity, probably due to the weak buffer capacity of shallow rooted grasslands and meadows to water deficit. For altitude, the vegetation vulnerability (AI) was higher at middle altitudes (2000–3000 m) than low altitudes (<2000 m) and high altitudes (3000–4500 m).

(3) During the period of 1982–2015, the AI was spatially heterogeneous. The 15 environmental factors that were investigated to evaluate the variation of AI accounted for 40% of the total variation based on the MLR. The relative contribution rate of silt soil fraction was the highest (23.3%). The water storage in areas with high silt soil content would be subjected to more evaporative losses, resulting in an increase of the vegetation vulnerability.