The TP is the largest alpine ecosystem in the world and has an extremely cold climate. The growing season temperature over the TP is about 12°C, which is similar to that in the Arctic (Huang et al., 2019). In contrast to the Arctic where the radiation is weak (~105 W m^-2) (Hamman et al., 2016), the TP receives strong solar radiation (~204 W m^-2), which is much stronger than that in the surrounding areas based on the observations from the meteorological stations of China Meteorological Administration (Yang et al., 2014). The precipitation regime over the TP is mainly shaped by the interaction between the Indian summer monsoon and the mid-latitude westerlies (Yao et al., 2012). The annual precipitation on the TP decreases from over 1000 mm in the southeast to less than 50 mm in the northwest (Kuang and Jiao, 2016), and the aridity status changes from humid and semi-humid in the southeast, to semi-arid in the central plateau, and arid in the northwest. Nearly 59.3% of the TP (1.5×10⁶ km²) is covered by alpine grasslands (Zhang et al., 2014). The grassland ecosystems vary along the aridity gradient from alpine meadow (0.7×10⁶ km²) in the southeast, to alpine steppe (0.58×10⁶ km²) in the northwest (Zhang et al., 2014).

Our analysis was conducted based on the framework proposed by Park et al. (2019). This framework analyzes climatic constrains for plant growth based on the law of minimum. It relays on two fundamental principles: (1) vegetation growth would be seasonally consistent with radiation if the ecosystem is not constrained by other climatic resources (e.g. temperature as a proxy for thermal energy, and precipitation as a proxy for water availability). Under this condition, the timing of peak vegetation growth would be close to the peak radiation. (2) When the ecosystem is constrained by thermal or hydrological climatic resources, the timing of peak vegetation growth would be close to the peak timing of the most limiting factor. This framework has been proven to be reliable, and was applied to investigate the climatic constraints for the peak photosynthesis over the northern hemisphere (Park et al., 2019; Zhao et al., 2022). Over the Tibetan Plateau, the timing of peak radiation occurs much earlier than the temperature and precipitation peaks. Due to the high elevation, the alpine ecosystem over the TP receives much stronger solar radiation than other temperate regions (Chen et al., 2014b), and radiation is not likely to be the limiting climatic resource for vegetation growth. Therefore, in this study we mainly focus on thermal and water constraints for peak vegetation growth. We also define a scenario when peak vegetation growth occurs between the peak temperature and peak precipitation as a transition state between thermal- and water-limited ecosystem ( Figure 1 ).

We used the normalized difference vegetation index (NDVI) to detect the peak growth date. The NDVI was derived from the Nadir Bidirectional Reflectance Distribution Function-Adjusted Reflectance (NBAR) from the MODerate-resolution Imaging Spectroradiometer (MODIS) (MCD43C4 v061). This latest version of the MODIS product does not have the major defect of sensor degradation that was present in the early datasets (Hilker et al., 2014) and should, therefore, be able to provide a more accurate estimate of changes in peak growth over the TP. The MCD43C4 dataset provides the reflectance for seven bands at daily time scale. We calculated NDVI from the reflectance (ρ) in the red and NIR bands as (ρNIR - ρRed)/(ρNIR + ρRed). We only used data flagged as “good quality”. The derived NDVI dataset has a spatial resolution of 0.05° × 0.05°, and cover the period from 2001 to 2020. We interpolated it into a spatial resolution of 0.1° × 0.1° for further analysis. We only considered vegetated regions, and pixels with annual NDVI below 0.1 was excluded as bare ground.

To test the robustness of our analysis, we also included four alternative vegetation indices and gross primary production (GPP) products to detect the peak growth date, including the NDVI derived from MODIS surface reflectance products (MOD09A1, v061), the GOSIF GPP product (Li and Xiao, 2019), the CSIF product (Zhang et al., 2018), and the PML-V2 GPP product (Zhang et al., 2019). We calculated NDVI from the MOD09A1 product for the period of 2001–2020, with a spatial resolution of 500 m and 8-day time step. We only used data with good quality, and records marked as cloud, cloud shadow, and aerosols were removed from analysis. GOSIF-GPP is the fine-resolution dataset developed from the OCO-2 based Global SIF product (GOSIF) and the SIF-GPP linear relationship at a 0.05° spatial resolution and 8-day time step (Li and Xiao, 2019). This GPP dataset covers the period 2001–2020. CSIF is generated by training artificial neural networks with MODIS surface reflectance and OCO-2-based solar induced chlorophyll fluorescence, with a 4-day temporal resolution and spatial resolution of 0.05° × 0.05° (Zhang et al., 2018). In this study, we used clear-sky instantaneous CSIF for the period 2001–2019. The PML-V2 dataset estimates GPP and ET simultaneously at 500m resolution and daily time scale using a coupled diagnostic biophysical model (PML-V2) forced by MODIS derived leaf area index, albedo, emissivity, and GLDAS meteorological datasets (Zhang et al., 2019). It was calibrated using 26 eddy flux (EC) stations in China. All these datasets were resampled to a spatial resolution of 0.1° × 0.1° and interpolated to a daily time step using spline interpolation at each pixel. We further used the PML-V2 GPP data to investigate the influence of the shift of peak growth date on ecosystem productivity.

The climate variables, including near-surface air temperature, precipitation, and incident solar radiation were derived from the reanalysis product of Indian Monsoon Data Assimilation and Analysis reanalysis (IMDAA) project. The IMDAA product is a regional reanalysis product that assimilates a variety of conventional observations including satellite retrievals and synthetic data though the Met Office Unified Model (UM) using four-dimensional variational (4D-Var) data assimilation method (Mahmood et al., 2018). It provides consistent dataset of high-resolution fields (0.12° × 0.12°) from 1979–2020 (Rani et al., 2021). We integrated the hourly dataset to a daily time step and interpolated to a spatial resolution of 0.1° × 0.1°.

We also analyzed in situ soil moisture from nine sites ( Table 1 ) to investigate the influence of permafrost thaw on the water availability for plant growth during peak season. These sites are located along the Qingzang railway, and the data was obtained from Zhao et al. (2021). We used the soil moisture observation at top layer of 10cm depth. For sites that lack observations at 10cm depth, we used records from the adjacent layer of 5cm or 15cm.

To calculate the peak vegetation growth date, we applied the singular spectrum analysis (SSA) to the daily NDVI or GPP within each year at each pixel, so as to remove high frequency noise and retain the relatively low frequency signals such as the seasonal cycle (Keenan et al., 2014; Zhou et al., 2017; Park et al., 2019). The peak growth date is then calculated as the date on which the filtered NDVI value reaches its annual maximum. The same method is applied to daily climatic variables (precipitation, temperature, and surface downward short-wave radiation) to identify the peak dates.

The date of peak vegetation growth over the TP occurs at early August (day of the year, DOY 213) ( Figure 2A ). It is months later than the date of peak radiation (DOY 125), and is close to the date of peak temperature (DOY 208) and peak precipitation (DOY 199) ( Figure 2A ), suggesting that the date of peak vegetation growth occurs around the period when the thermal and hydrological condition is most suitable for vegetation growth. During this period, the snow cover over the TP reaches its minimum ( Figure 2A ) and the top-layer soil is thawed due to the warm climate in summer (Wu and Zhang, 2010). This multi-year averaged peak growth date suggests a long-term acclimation of vegetation growth to the regional climate over the TP, and it supports our assumption that the date of peak growth occurs around the peak of most constraining factor.

Spatially, the date of peak vegetation growth delays from late-July in the east (DOY 215 ± 5) to mid-August (DOY 228 ± 2) in the west along the aridity gradient ( Figure 2B ). In the humid and semi-humid region where the annual precipitation is reaches 250 mm year^-1, the peak precipitation (DOY 197) occurs 11 ± 5 days before the peak temperature (DOY 208), and the vegetation growth reaches the peak at 10 ± 7 days later than the date of peak temperature, suggesting that the alpine meadow in the humid and semi-humid region is constrained by temperature. In the semi-arid region with an annual precipitation below 200 mm year^-1, the peak precipitation lags behind peak temperature for 8 ± 3 days, and the vegetation growth peaks at 20 ± 2 days later than the peak precipitation. This indicates that the water availability is the primary constrain for the ecosystem in this semi-arid region. Spatially, 46% of the region is thermal-constrained, with 77% of which located in the semi-humid east, and 48% of the region is water-constrained, with 85% of which lay in the semi-arid west. For 6% of the region, the timing of peak vegetation growth occurs between the peak temperature and peak precipitation, suggesting it in a transition state of thermal- and water-constrained ecosystem ( Figure 2C ). This spatial pattern of peak growth date and the relative positioning against the climatic peaks demonstrated that our analysis framework well represents the primary constrains of vegetation growth over the plateau.

The date of peak vegetation growth over the TP shifted notably during 2001–2020, while the shifting patterns are divergent between the eastern and western plateau ( Figure 3A ). In the east, 58% of the grasslands showed significant advancing peak growth. On average, the peak growth date of eastern plateau significantly advanced at a rate of 3.72 ± 1.36 days decade^-1, which keeps in pace with the advancing date of the onset of growing season (Shen et al., 2022). In the west, however, the peak vegetation growth delayed at a rate of 0.37 ± 0.64 days decade^-1, with 20% of the region showed a significant delaying trend. On average, the date of peak vegetation growth over the TP advances at a rate of 0.81 ± 0.71 days decade^-1 during 2001–2020 ( Figure 3B ). Results from other four vegetation indices and GPP products showed similar shifting patterns, demonstrated the robustness of the changes in peak growth dates over the TP ( Supplementary Figure 1 ).

We then analyzed the changes of the relative positioning between the date of peak vegetation growth and that of climate peaks. The relative timing between the dates of peak vegetation growth and peak temperature (δ_G,T ) has changed significantly over the past 20 years, with significant increasing trend in the east and decreasing trend in the west, while that against the date of peak precipitation (δ_G,Pr ) is much less notable ( Figure 4 ). We then specifically analyzed the relative positioning of peak vegetation growth and climate peaks for 2000s (2001–2010) and 2010s (2011–2020) for the eastern and western plateau respectively ( Figures 5A–D ). In the eastern plateau, the climate shifts as an “advancing peak precipitation – delaying peak temperature” pattern. The date of peak precipitation advanced by 6 days from 198 ± 12 in 2000s to 192 ± 5 in 2010s, and the date of peak temperature varies delays from 207 ± 5 in 2000s to 210 ± 8 in 2010s. In response to these shifts in climate, the date of peak vegetation growth tracks the advancing peak precipitation and shifted spring-ward for 6 days from 219 ± 6 in 2000s to 213 ± 7 in 2010s ( Figures 5B, D ). In the western plateau, however, the climate shifts as a pattern of “stable peak precipitation – advancing peak temperature”. The date of peak precipitation changed by 1 day from 207 ± 7 in 2000s to 208 ± 9 in 2010s, and the date of peak temperature advanced by 5 days from 203 ± 6 in 2000s to 198 ± 11 in 2010s. The date of peak growth changes slightly from 225 ± 6 in 2000s to 226 ± 7 in 2010s ( Figures 5A, C ). As a consequence, δ_G,Pr decreased slightly by 1 day, while δ_G,T had a large increase of 8 days between 2000s and 2010s. To sum up, 10% of the grasslands over the TP turned from thermal-constrained ecosystem in 2000s into transition state from thermal- to water- constrained in 2010s, and most of these transitions occur in the semi-humid east ( Figures 5E, F ). Similar transition was also found when using different vegetation indices and GPP products to detect peak growth date ( Supplementary Figure 2 ). These shifts in the relative positioning between peak vegetation growth and climate peaks suggest that the alpine grassland ecosystem over the TP is moving towards the tipping point of a transition from thermal- to water-constrained ecosystem.

We found a notable “earlier peak – larger production” changes for alpine grasslands over the TP. In the semi-humid eastern TP, the ecosystem production increases by 1.99 gCm^-2 year^-1 per one day advances in peak vegetation growth, and this sensitivity decreases to 0.29 gCm^-2 year^-1day^-1 in the semi-arid west. The influence of the shifting peak vegetation growth on ecosystem production also changes over time. In the eastern plateau, the sensitivity of ecosystem production to shifting peak vegetation growth (β) decreases from -3.02 gCm^-2 year^-1day^-1 during 2000s to -1.25 gCm^-2 year^-1day^-1 during 2010s. In the western plateau, the β varies from -0.13 gCm^-2 year^-1day^-1 to -0.62 gCm^-2 year^-1day^-1 for the last two decades and is not significant for both periods ( Figure 6 ). This result is robust to the used of GPP dataset ( Supplementary Figure 3 ).