Heatwaves Signiﬁcantly Slow the Vegetation Growth Rate on the Tibetan Plateau

: 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 signiﬁcant trends from 2000 to 2020 for the entire Tibetan Plateau. On a monthly scale, the average frequency of heatwaves increased signiﬁcantly ( p < 0.1) in August, while no signiﬁcant 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.


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][5][6][7][8][9][10][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][15][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 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.

Study Area
The Tibetan Plateau, with an area of 2.5 million km 2 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.

Datasets
The daily maximum temperature (T max ) 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.

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: where G = 2.5, C 1 = 6, C 2 = 7.5, and L = 1; ρ blue , ρ red , and ρ nir are the reflectance of the blue, red, and near-infrared bands, respectively.
where Y represents the original time-series data, Y* is the reconstructed time-series data, C i 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].
where VIs bf and VIs dur are the average of VIs before the heatwave and during the heatwave, respectively. VGR baseline represents the multi-year average VGR during the reference period (2000-2020). ∆VGR is the difference between VGR and VGR baseline . 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.

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.  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.

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. 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. 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. 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).

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][49][50][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. 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.

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.
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.