*Article* **Heatwaves Significantly Slow the Vegetation Growth Rate on the Tibetan Plateau**

**Caixia Dong 1,2, Xufeng Wang 1,\*, Youhua Ran <sup>1</sup> and Zain Nawaz <sup>1</sup>**


**\*** Correspondence: wangxufeng@lzb.ac.cn; Tel.: +86-931-4967724

**Abstract:** In recent years, heatwaves have been reported frequently by literature and the media on the Tibetan Plateau. However, it is unclear how alpine vegetation responds to the heatwaves on the Tibetan Plateau. This study aimed to identify the heatwaves using long-term meteorological data and examine the impact of heatwaves on vegetation growth rate with remote sensing data. The results indicated that heatwaves frequently occur in June, July, and August on the Tibetan Plateau. The average frequency of heatwaves had no statistically significant trends from 2000 to 2020 for the entire Tibetan Plateau. On a monthly scale, the average frequency of heatwaves increased significantly (*p* < 0.1) in August, while no significant trends were in June and July. The intensity of heatwaves indicated a negative correlation with the vegetation growth rate anomaly (ΔVGR) calculated from the normalized difference vegetation index (NDVI) (r = −0.74, *p* < 0.05) and the enhanced vegetation index (EVI) (r = −0.61, *p* < 0.1) on the Tibetan Plateau, respectively. Both NDVI and EVI consistently demonstrate that the heatwaves slow the vegetation growth rate. This study outlines the importance of heatwaves to vegetation growth to enrich our understanding of alpine vegetation response to increasing extreme weather events under the background of climate change.

**Keywords:** heatwave; alpine vegetation; Tibetan Plateau; remote sensing; extreme climate events

#### **1. Introduction**

A heatwave is defined as a period with sustained high-temperature anomalies resulting in strong impacts on human health, the ecological environment, and socioeconomic development [1]. A recent study has indicated that heatwaves have increased in prevalence significantly since the 1950s [2]. The heatwave has received growing attention in global change ecology study because of its remarkable effects on carbon, water, and energy exchange between the land surface and atmosphere [3]. Much evidence from crop yield, tree ring, and manipulative experiments has demonstrated that the occurrence of extreme hightemperature events can trigger significant impacts on terrestrial ecosystems and human society [4–11]. The heatwave also significantly impacts vegetation from regional to global scales, which has been witnessed by the satellite-derived normalized difference vegetation index (NDVI)-based studies [12,13]. These studies have mainly focused on tropical and temperate regions but have ignored cold regions.

With global warming, cold regions have been experiencing increased intensity, frequency, and duration of heatwaves in the past several decades [14–16]. As the "Third Pole", the Tibetan Plateau is warming twice as fast as the global average warming rate [17]. In situ meteorological observations and model projections indicated that extreme hightemperature events have happened frequently on the Tibetan Plateau. Due to the low intensity of anthropogenic activities [18], the Tibetan Plateau is an ideal region for studying the responses of alpine vegetation to extreme temperature events. The alpine vegetation on the Tibetan Plateau is very sensitive to high-temperature events due to the heat-limiting

**Citation:** Dong, C.; Wang, X.; Ran, Y.; Nawaz, Z. Heatwaves Significantly Slow the Vegetation Growth Rate on the Tibetan Plateau. *Remote Sens.* **2022**, *14*, 2402. https://doi.org/ 10.3390/rs14102402

Academic Editor: Miaogen Shen

Received: 8 April 2022 Accepted: 13 May 2022 Published: 17 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

environment [19]. Many recent studies have begun to examine the impacts of extreme temperature on vegetation (e.g., productivity, phenology) on the Tibetan Plateau [20,21]. As a specific type of extreme high–temperature event, the heatwave has been reported recently by literature [22] on the Tibetan Plateau. However, very few studies have been conducted on the ecological effects of the heatwave on the Tibetan Plateau. We only found one site scale study, which reported that the heatwave can substantially increase alpine ecosystem respiration on the Tibetan Plateau [20]. Therefore, the response of alpine vegetation on the Tibetan Plateau to heatwaves is poorly understood. It is necessary to evaluate the heatwave effect on vegetation more widely and provide valuable information to address the climate change in this region.

Meanwhile, the MODIS bidirectional reflectance distribution function (BRDF)-adjusted daily reflectance product has made it possible to detect vegetation change in a very short period. On the Tibetan Plateau, heatwaves are usually of short duration and take place at a small part of the plateau. Some small heatwaves may be ignored if using the 16-day or monthly composite vegetation index. The BRDR-corrected daily snow–free MODIS reflectance product and long-term daily meteorological observations make it possible to examine the alpine vegetation response to heatwaves at a regional scale on the Tibetan Plateau.

We aim to fill the knowledge gap about the response mechanism of the alpine vegetation to heatwave on the Tibetan Plateau. The objectives of this study include: (1) identify the heatwaves from 2000 to 2020 on the Tibetan Plateau and analyze their temporal trends; (2) examine the response of vegetation growth to heatwaves intensity and duration.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Tibetan Plateau, with an area of 2.5 million km2 and an average altitude over 4000 m above sea level [23], is the largest and highest plateau in the world [4]. The geographical range of the Tibetan Plateau is 26–40◦N, 73–105◦E. It spans six provinces, namely the Tibet and Xinjiang Uygur Autonomous Regions and Qinghai, Yunnan, Sichuan, and Gansu Provinces [15]. The characteristics of the climate on the Tibetan Plateau are strong solar radiation, low air temperature, and large day–night temperature difference [18,24]. This climatic pattern determines the general distribution of the vegetation [25]. The dominant vegetation type is alpine grasslands, which accounts for about 60% of the entire plateau area [26]. Figure 1 depicts the spatial distribution of the meteorological stations on the Tibetan Plateau.

**Figure 1.** The distribution of the meteorological stations on the Tibetan Plateau.

#### *2.2. Datasets*

The daily maximum temperature (Tmax) and precipitation data used in this study were provided by the Climatic Data Center, National Meteorological Information Center, China Meteorological Administration (http://data.cma.cn/, accessed on 21 August 2021). In the dataset, there are 86 meteorological stations across the Tibetan Plateau. Because heatwave detection requires long-term temperature data, we excluded the sites with data gaps from 1980 to 2020. Finally, 64 stations were selected in this study and are depicted in Figure 1. Table A1 in Appendix A depicts information about each meteorological station (WMO code, name, latitude, longitude, and elevation.

Digital elevation model data for the Tibetan Plateau was obtained from Shuttle Radar Topography Mission (https://earthexplorer.usgs.gov/, accessed on 21 August 2021) and its spatial resolution is 30 m.

The MODIS Nadir Bidirectional Reflectance Distribution Function Adjusted Reflectance (NBAR) product (MCD43A4) was used to monitor the vegetation growth in this work. MCD43A4 provides daily surface reflectance by combining Terra and Aqua MODIS data at a 500–m spatial resolution. The surface reflectance was normalized to nadir using the bidirectional reflectance distribution function for the solar angle at local noontime. This product has removed view angle effects, and masks cloud cover and snow contamination [27]. Collection 6 of MCD43A4 from 2000 to 2020 on the Tibetan Plateau was obtained from the Google Earth Engine (https://earthengine.google.com/, accessed on 21 August 2021). The two vegetation indices (Vis) have been widely used in ecological studies, whereas the NDVI is chlorophyll-sensitive, the EVI is more responsive to canopy structural variations [28,29]. A combination of the NDVI and EVI to complement each other examined the robustness and comparability of the vegetation growth rate changes [30]. In this study, we used the surface reflectance from the MCD43A4 product to calculate daily Vis, including the NDVI and EVI.

The Land Surface Soil Moisture Dataset over the Tibetan Plateau was downloaded from National Tibetan Plateau Data Center [31]. This dataset is daily surface soil moisture with a spatial resolution of 0.25◦, retrieved from passive microwave brightness temperature data. The dataset synthesized microwave brightness temperature measurements from SMMR, SSM/I, SSMIS, AMSR-E, AMSR2, SMAP, and FY3B to produce a long-term soil moisture product [32]. We used Land Surface Soil Moisture Dataset data from 2002 to 2020 to examine the soil moisture difference before and during the heatwave.

#### *2.3. Methods*

In this study, we used a percentile-based thresholds method to identify the heatwaves in June, July, and August at each station on the Tibetan Plateau. [33]. The heatwaves were identified following the definition of a heatwave in the literature [34] with a few modifications in this study. The 90th percentile values of daily maximum temperature during the climatological baseline period (1980–2020) are used as the threshold to identify hot days. A period with at least five consecutive hot days (the maximum temperature is greater than the threshold) was identified as a heatwave event at a single station. The frequency, duration, and intensity are characteristics of a heatwave event. The frequency of heatwaves is the sum of the heatwaves at 64 sites on the Tibetan Plateau in a year. The number of stations experiencing heatwaves is also counted each year. Each site is counted only once per year. The duration and temperature anomaly above the threshold are introduced as two dimensions to quantify the severity of the heatwaves. The duration is the length of a heatwave event, and the temperature anomaly is the average temperature difference between the daily maximum temperature and the 90th percentile threshold. The accumulative intensity of a heatwave is the sum of temperature differences above the threshold during the heatwave. The average intensity of heatwaves is the average temperature difference during the heatwave. To explore the impact of heatwaves on vegetation, we only focused on the heatwaves occurring in June, July, and August.

VIs were used to calculate the vegetation growth rate, which is the change of the value of VIs before and after a heatwave. The NDVI and EVI [30] are calculated with the Equations (1) and (2), respectively. In addition to correct inferior values in VIs, a time series reconstruction algorithm, the Savitzky–Golay filter (Equation (3)), is applied to long-term daily VIs in this study [35]. We excluded these heatwaves in the analysis when MCD43A4 data was missing to make our result reliable. Vegetation growth rate (VGR) was calculated as the difference between the vegetation index before and during the heatwave (Equation (4)). Then, the vegetation growth rate anomaly (ΔVGR) (Equation (5)) was calculated from the VGR during the heatwave and the multi–year average VGR in the corresponding period. The vegetation index and vegetation growth rate anomaly were used as vegetation proxies to examine the heatwave impacts on vegetation. The formulas are expressed as follows:

$$\text{NDVI} = \frac{\rho\_{nir} - \rho\_{red}}{\rho\_{nir} + \rho\_{red}} \tag{1}$$

$$\text{EVI} = \frac{G \times (\rho\_{nir} - \rho\_{rel})}{\rho\_{nir} + (\mathbf{C}\_1 \times \rho\_{nir} - \mathbf{C}\_2 \times \rho\_{bluc}) + L} \tag{2}$$

where *G* = 2.5, *C*<sup>1</sup> = 6, *C*<sup>2</sup> = 7.5, and *L* = 1; *ρblue*, *ρred*, and *ρnir* are the reflectance of the blue, red, and near-infrared bands, respectively.

$$Y\_j^\* = \left(2m+1\right)^{-1} \sum\_{i=-m}^{i=m} C\_i Y\_{j+i} \tag{3}$$

where *Y* represents the original time–series data, *Y\** is the reconstructed time-series data, *Ci* is the weight of the filter window, and *2m* + 1 is the size of the filter window. The window size and polynomial order in the Savitzky–Golay filter were set to 31 and 4, respectively [36].

$$\text{VGR} = VIs\_{dur} - VIs\_{bf} \tag{4}$$

$$
\Delta \text{VGR} = VGR - VGR\_{\text{baseline}} \tag{5}
$$

where *VIsbf* and *VIsdur* are the average of *VIs* before the heatwave and during the heatwave, respectively. *VGRbaseline* represents the multi-year average VGR during the reference period (2000–2020). ΔVGR is the difference between VGR and *VGRbaseline*.

The linear trend of the number of sites with heatwaves was analyzed using the Mann– Kendall methods [37,38]. The Mann–Kendall method is a nonparametric test for monotonic trends. This method does not assume a specific distribution for the data and is not sensitive to outliers. The Theil–Sen method was used to calculate the slope of the linear trend [39]. The slope of the trend measures the number of the heatwaves' change rate over time. To explore the impact of heatwaves on vegetation growth, we calculated the correlation coefficients between the heatwaves (the intensity of the heatwaves) and vegetation growth rate (change rate anomaly of NDVI and EVI) using the Pearson correlation method.

#### **3. Results**

#### *3.1. Trends of Heatwaves Frequency*

Based on the 64 meteorological stations on the Tibetan Plateau, we first identified the heatwaves and calculated the duration and intensity of the heatwaves for each station. The interannual variation of heatwaves frequency at these stations is depicted in Figure 2. From 2000 to 2020, the heatwaves frequently occurred in June, July, and August. Overall, from 2000 to 2020, the frequency of heatwaves in the growing season had no statistically significant trends for the entire Tibetan Plateau in Figure 2. The occurring frequencies of the heatwaves are different among June, July, and August. The heatwaves happened more frequently in August than in June and July. The heatwave frequency increased significantly (*p* < 0.1) in August, while no significant trends occurred in June and July.

**Figure 2.** The average frequency of heatwaves at 64 stations in June (**a**), July (**b**), and August (**c**), and June−August (**d**) on the Tibetan Plateau from 2000 to 2020.

Figure 3 depicts the interannual dynamic of the duration and intensity of heatwaves from 2000 to 2020. The color changes from yellow to red indicate that the heatwave severity varies from weak to strong. Overall, the duration and intensity of heatwaves ranged from 5.00 to 11.57 days and from 0.67 to 2.55 ◦C/d, respectively. Most of the heatwaves are short with a duration of about 6 days. The longest duration appeared in August. Among June, July, and August, heatwaves in August lasted longer than in other months. The intensity of the heatwaves has neither obvious monthly patterns nor evident temporal changing trends. Through analyzing the intensity and duration of heatwaves, we found that heatwaves with long durations may have low intensity (average high–temperature anomaly). The heatwaves occurred most frequently in recent years, such as 2006, 2013, 2016, etc.

To explore the extent of the heatwaves, a matrix heatmap is used to depict the number of stations where heatwaves happened in June, July, and August from 2000 to 2020 (Figure 4). Generally, as expected, there are no widespread heatwave occurrences in most years on the Tibetan Plateau due to the cold environment. From 2000 to 2020, heatwaves occurred at more than half of the total stations on the Tibetan Plateau in 2000, 2001, 2006, 2010, 2013, and 2016. The most two severe heatwave events detected from the daily maximum temperature happened in August 2016 with 47 sites and in June 2013 with 38 sites.

**Figure 3.** Matrix heatmap for heatwave duration and intensity in June, July, and August on the Tibetan Plateau from 2000 to 2020. The matrix heatmap (**a**) refers to the temporal change of the heatwave duration; the matrix heatmap (**b**) refers to the temporal change of the heatwave intensity. Each grid cell represents the average duration and intensity of a heatwave at all sites in a given month of a year. The blank grid cell represents no heatwave, and the value is NaN. The color of the matrix heatmap represents the size of the value.

**Figure 4.** Matrix heatmap of the number of sites with heatwaves from June to August on the Tibetan Plateau from 2000 to 2020. Each grid cell represents the number of sites with heatwaves in a given month of a year. The blank grid represents no heatwave, and the value is NaN. The color of the matrix heatmap represents the size of the value.

#### *3.2. Effects of Heatwaves on Vegetation*

The deviation analysis was used to calculate the ΔVGR, which can reflect the VGR change caused by the heatwave. A positive anomaly indicates an increase in the VGR, and a negative anomaly indicates a decrease in VGR. Figure 5 depicts the average ΔVGR calculated from NDVI (a) and EVI (b) in June, July, and August from 2000 to 2020. The ΔVGR ranging from positive to negative was expressed as the color changing from green to yellow. The range of the ΔVGR calculated by NDVI and EVI is from −0.0088 to 0.057 and from −0.0095 to 0.023, respectively. Overall, the ΔVGR calculated from NDVI is consistent with the ΔVGR calculated from EVI.

**Figure 5.** Matrix heatmap of the ΔVGR from June to August on the Tibetan Plateau from 2000 to 2020. The matrix heatmap (**a**) refers to the ΔVGR calculated from NDVI; the matrix heatmap (**b**) refers to the ΔVGR calculated from EVI. Each grid cell represents the average vegetation growth rate anomaly at all sites in a given month of a year. The blank grid cell represents no heatwave happened, and the value is NaN. The color of the matrix heatmap represents the size of the value.

The monthly variation of ΔVGR corresponding to heatwaves on the Tibetan Plateau in June, July, and August is displayed in Figure 6. The ΔVGR was negative when the heatwave happened in June, July, and August. That means the VGR during the heatwaves is lower than the multi-year mean on the Tibetan Plateau. Compared to other months, the minimum value of the ΔVGR was found in June, indicating that vegetation growth in June was more sensitive to the heatwave than in July and August. The value of the ΔVGR calculated by NDVI and EVI corresponding to the heatwaves is more consistent in June and July than in August. Overall, heatwaves in the growing season significantly slow the VGR on the Tibetan Plateau. The analysis indicates that the ΔVGR can capture the vegetation growth response to heatwaves on the Tibetan Plateau.

**Figure 6.** The ΔVGR in each month in June, July, and August on the Tibetan Plateau from 2000 to 2020. The ΔVGR is calculated from NDVI (**a**) and EVI (**b**) on the Tibetan Plateau from June to August. The black lines represent the monthly multiyear mean of each variable between 2000 and 2020; the red lines indicate the seasonal evolution during extremely high temperatures.

The correlation relationships were analyzed between the ΔVGR and the intensity of the heatwaves. The intensity of the heatwaves was calculated as temperature anomaly during heatwaves multiplied by heatwave duration. The intensity was grouped with a step of 2.5 ◦C×d to make the analysis clearer. As depicted in Figure 7a,c, the accumulative intensity (◦C×d) of heatwaves indicates a negative correlation with the ΔVGR calculated from NDVI (r = −0.74, *p* < 0.05) and EVI (r = −0.61, *p* < 0.1) on the Tibetan Plateau, respectively. The average intensity (◦C/d) of heatwaves indicates a negative correlation with the ΔVGR calculated by NDVI (r = −0.77, *p* < 0.05) and EVI (r = −0.66, *p* < 0.1) on the Tibetan Plateau (Figure 7b,d), respectively. Vegetation growth is strongly affected by the heatwaves in June, July, and August on the Tibetan Plateau. With the intensity increase of heatwaves, the VGR decreases linearly.

**Figure 7.** Pearson's correlation coefficient between the intensity of the heatwave and the ΔVGR from June to August during the growing season for all 64 sites on the Tibetan Plateau. (**a**,**c**) the accumulative intensity of heatwaves versus the ΔVGR. (**b**,**d**) the average intensity of heatwaves versus the ΔVGR. The ΔVGR in (**a**,**b**) is calculated from NDVI, and in (**c**,**d**) is calculated from EVI.

To corroborate our findings, we select the year 2013 and 2016, when widespread heatwaves occurred, to specially study the anomaly of VGR and the anomaly of climate factors before and after the occurrence of the heatwaves. The spatial distribution of the ΔVGR on the Tibetan Plateau is displayed in Figure 8a for June 2013 and Figure 8c for August 2016. The average VGRs in June 2013 and August 2016 were significantly lower than the multi-year average VGR. In June 2013, the ΔVGRs were negative at 33 of the 38 sites where the heatwave occurred. In August 2016, the ΔVGRs were negative at 30 of the 47 sites where the heatwave occurred. An obvious decrease in vegetation growth rate can be found in most sites where the heatwave occurred. It is suggested that the heatwaves in 2013 and 2016 significantly slowed down the VGR. During the selected two heatwave events, the temperature was significantly higher than the multi−year average in the corresponding period, including 5.8 ◦C above the multi−year average in 2016 and 4.9 ◦C above the multi−year average in 2013 (Figure 9).

**Figure 8.** ΔVGR response to heatwaves in June 2013 (**a**) and August 2016 (**c**). The blank circle represents the site without a heatwave, while the filled circle represents the site with a heatwave; (**b**,**d**) represent the site number in June 2013 and August 2016, respectively; the red-filled circle represents the negative ΔVGR; the green−filled circle represents the positive ΔVGR. Positive ΔVGR stands for an increase in the VGR; negative ΔVGR stands for a decrease in the VGR. The right plots illustrate the number of sites with negative and positive VGR.

**Figure 9.** The ΔVGR and daily maximum temperatures during the heatwaves in June 2013 and August 2016. (**a**,**c**) represent the ΔVGR in June 2013 and August 2016, respectively; (**b**,**d**) represent the average daily maximum temperature (Tmax) versus the multi-year average Tmax during the heatwave occurrence in June 2013 (**b**) and August 2016 (**d**). The red line indicates the Tmax (15 June 2013– 20 June 2013 DOY: 166–171 and 16 August 2016–24 August 2016 DOY: 228–236) during a heatwave in one severe year. The black line indicates the multi-year average Tmax during heatwaves occurrence in June 2013 and August 2016.

#### **4. Discussion**

To our knowledge, very few studies have focused on heatwaves on the Tibetan Plateau. Previous studies mainly used the traditional extreme high-temperature indices to explore their effects on alpine vegetation, such as TX90p (percentage of days when TX > 90th percentile), WSDI (warm spell duration index), etc. The traditional extreme high-temperature indices are usually calculated at a monthly or annual scale, which is too coarse to capture the short climate disturbance on vegetation. A recent study examined the trends of extreme temperature events using 71 meteorological stations from 1961 to 2005 and found that there were statistically significant increasing trends for extreme high-temperature indices [15]. He et al. [21] analyzed the spatial pattern and long-term trend in extreme high-temperature indices in the Kobresia meadow region from 1961 to 2008, and found a significant increase in the warmest daytime temperature. However, in this study, we found that heatwaves have no significant increasing trend from 2000 to 2020, and the trends vary greatly among June, July, and August. This is inconsistent with the extreme high−temperature indices trend reported previously on the Tibetan Plateau [40]. This is caused by the different definitions between heatwave and extreme high-temperature events. The inconsistency in trends between heatwave and extreme high-temperature indices could also be attributed to the different periods and datasets. Liu et al. [40] reported the characteristic of extreme high-temperature events from 2001 to 2015 based on the monthly gridded datasets, but this study explored the heatwaves based on the meteorological station data. The heatwaves on the Tibetan Plateau mainly occurred locally, only a few widespread heatwaves were detected, and no heatwaves occurred for the entire Plateau or all sites simultaneously. Due to the sparse and non-uniformly distributed weather stations, it is difficult to accurately

extract the heatwave spatial extent [41]. However, the evolution of spatial extent is essential to better understand the varying mechanism of heatwaves on the Tibetan Plateau. Thus, future studies on better understanding the dynamic of heatwaves will be benefited from the high-resolution and reliable grid meteorological dataset.

Heatwaves can limit vegetation photosynthesis by pushing the ambient temperature to exceed the optimal photosynthetic temperature, increasing the vapor pressure deficit (VPD), and reducing the soil moisture. High temperatures over the optimal photosynthetic temperature could constrain the activity of Rubisco, increase photorespiration, and lead to a decline in net photosynthesis [42]. An ambient temperature lower or higher than the optimal photosynthetic temperature will inhibit vegetation growth [43]. The slowdown of vegetation growth rate during heatwave occurrence indicates that the extreme temperatures of heatwaves overpassed the optimum photosynthesis temperature for alpine vegetation on the Tibetan Plateau. Additionally, summer heatwaves affect photosynthesis primarily due to the physiological response to water deficit and high temperatures, including reductions in enzymatic activity and stomatal conductance to prevent water loss [44]. The stress effects are increased by water deficits [45]. In general, the occurrence of heatwaves was frequently accompanied by a decline in precipitation and a decrease in air relative humidity [5]. To validate this phenomenon, we examined the soil moisture and precipitation for the two selected heatwaves in June 2013 and August 2016. The soil moisture and precipitation during heatwave periods are obviously lower than the multi−year average in the corresponding periods (Figure 10). This demonstrates that heatwaves affect alpine vegetation by combining temperature stress and water limitation on the Tibetan Plateau. However, comparing the soil moisture and precipitation before and during the heatwave in August 2016, soil moisture and precipitation showed different change patterns. Maybe, this resulted from data noises in soil moisture product. The microwave–based soil moisture used in this work is retrieved using in–situ at five pixel-scale fields and 25km microwave remote sensing [32]. But, it is not widely validated on the Tibetan Plateau due to no widespread in-situ measurements. Moreover, the alpine vegetation responds differently to heatwaves in different phenology stages. The ΔVGR in July is more significant than that of July and August. It is indicated that the alpine vegetation is more sensitive to heatwaves in the early growing season than in the later growing season. Vegetation is fragile and sensitive to the environment in the early growing season [26,45,46], for example, spring phenology is more sensitive to environmental factors than autumn phenology [47]. Meanwhile, vegetation usually grows faster in the early growing season than in the later growing season; therefore, the growth rate may be more sensitive to environmental factors in the early growing season [46]. In this study, it is indicated that ΔVGR decreased linearly with heatwave intensity. This is partly because most heatwaves are weak on the Tibetan Plateau, and alpine vegetation can recover from these disturbances. Different grasslands on the Tibetan Plateau exhibited different response patterns to climate changes [26]. Interestingly, the alpine meadow is more sensitive to heatwaves than the alpine steppe in June (Figure 11), but the opposite is true in July and August. This may result from the different coverage between the two types. Moreover, the growth rate changing mechanism is complex; more factors should be considered [48–51]. Further research is needed to clarify the detailed mechanism of these changes. On the Tibetan Plateau, the vegetation that has adapted to the cold, alpine environment might differ from other ecosystems concerning vegetation responses to temperature extremes [26,40]. Given that the Tibetan Plateau will continue warming in the future [52,53], heatwaves will happen more widely and intensely and will lead to abrupt changes by approaching the temperature threshold of alpine vegetation.

**Figure 10.** Precipitation and soil moisture comparison before and during the heatwave occurrence in June 2013 and August 2016. (**a**,**c**) represent the average soil moisture before/during two selected heatwave periods versus multi-year average soil moisture in the corresponding periods. (**b**,**d**) represent the average daily precipitation before/during two selected heatwave periods versus multi-year average soil moisture in the corresponding periods. The red bars indicate the heatwave years. The black bars indicate the baseline period (soil moisture: 2002–2019, pre: 2000–2020). The "bf" represents the period before the heatwaves (9 June 2013–14 June 2013 DOY: 160–165 and 9 August 2016–15 August 2016 DOY: 221–227). The "dur" indicates the heatwave periods (15 June 2013–20 June 2013 DOY: 166–170 and 16 August 2016–24 August 2016DOY: 228–236).

There are some uncertainties involved in this study. Firstly, meteorological observations have a relatively sparse and uneven distribution, resulting in the low representativeness of the identification of heatwaves [41]. Secondly, the quality of the remote sensing dataset is low on the Tibetan Plateau due to the contamination of snow, clouds, and complex terrain. Thirdly, the reconstructed vegetation index value can result in uncertainties in the analysis. The different spatial representativeness among station data, MODIS data, and coarse resolution soil moisture data can also lead to uncertainties in the study. Fourthly, the definition of a heatwave is uncertain due to the unique natural conditions on the Tibetan Plateau. In some heatwave definitions in the tropic or temperate region, a fixed high-temperature threshold is usually used by combining the temperature 90% percentiles. In this work, only the 90% percentile was used and may result in uncertainties when comparing heatwaves on the Tibetan Plateau with other regions.

**Figure 11.** ΔVGR in each month in June, July, and August on the Tibetan Plateau from 2000 to 2020. The dark green column represents the ΔVGR of the alpine meadow; the light green column indicates the ΔVGR of the alpine steppe.

#### **5. Conclusions**

In this study, heatwaves were detected on the Tibetan Plateau by using long-term meteorological station data. The characteristics of heatwaves were explored at these meteorological stations. By combining the remotely sensed vegetation indices, the alpine vegetation response mechanism to heatwaves was examined on the Tibetan Plateau. The results indicate that: (1) With rapid warming, heatwaves frequently occur in June, July, and August on the Tibetan Plateau. The heatwaves have no significant increasing trend from 2000 to 2020; (2) The correlation between heatwave intensity and vegetation growth rate anomalies was significantly negative on the Tibetan Plateau. The vegetation growth rate estimated from NDVI and EVI consistently indicates that heatwaves slow vegetation growth. The alpine vegetation growth rate is more sensitive to the heatwave in June than in July and August. This study outlines the importance of heatwaves to vegetation growth to enrich our understanding of alpine vegetation response to increasing extreme weather events under the background of climate change.

**Author Contributions:** C.D. and X.W. Data collection, conceptualization, writing—original draft preparation, software, methodology, formal analysis; X.W. funding acquisition, writing review and editing, supervision; Y.R. and Z.N. writing—reviewing and editing the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (Grant No. 41771466), the Youth Innovation Promotion Association CAS to X.W. (NO. 2020422), and the National Key R&D Program of China (Grant No. 2017YFA0604801).

**Data Availability Statement:** The meteorological data is available at the Climatic Data Center, National Meteorological Information Center, China Meteorological Administration (http://data. cma.cn/, accessed on 7 April 2022). The DEM (Digital Elevation Model) data is available at http: //srtm.csi.cgiar.org/, accessed on 7 April 2022. The Collection 6 of MCD43A4 is available in the Google Earth Engine (https://earthengine.google.com/, accessed on 7 April 2022). The soil moisture data is available at the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn/, accessed on 7 April 2022).

**Acknowledgments:** The authors would like to thank TPDC for providing the data and anonymous reviewers for their valuable comments.

**Conflicts of Interest:** The authors declare no competing interest.

#### **Appendix A**

**Table A1.** Information for meteorological stations used in this study.



**Table A1.** *Cont.*

#### **References**


#### *Article*
