Duration of vegetation green-up response to snowmelt on the Tibetan Plateau

The Tibetan Plateau (TP) is characterized by abundant snow and heightened sensitivity to climate change. Although the impact of snowmelt on vegetation green-up is well recognized, the duration of vegetation response to snowmelt on the TP remains unclear. This study calculates the time differences between the green-up date and the start of snowmelt from 2001–2018 on the TP, denoted as ΔD. Exploratory spatial data analysis and the Mann–Kendall test were then applied to investigate the spatiotemporal distribution feature of ΔD. Subsequently, partial correlation and multiple linear regression analysis were employed to examine the impact of spring mean temperature, spring total rainfall, and daily snowmelt on ΔD. The results reveal that the mean ΔD across the TP was 38.5 d, with a spatially clustered distribution: low ΔD values were concentrated in the Hengduan Mountains, while high ΔD values were observed in the Bayankara and Himalaya mountains. Additionally, ΔD shortened with increasing spring temperature, total rainfall, and daily snowmelt, which accounted for 23.5 %, 28.8 %, and 35.4 % of ΔD variation, respectively. In 67 % of arid areas and 64 % of regions with low vegetation, daily snowmelt was the dominant factor influencing ΔD. Conversely, spring temperature played a primary role in 48 % of humid areas and 37 % of regions with high vegetation. Our findings enhance the understanding of vegetation responses to snowmelt and provide a scientific foundation for further research on the stability of alpine ecosystems and the impacts of climate change on the TP.

1 Introduction

The Tibetan Plateau (TP) is the largest plateau in China and the highest in the world. It serves as the source of several major Asian rivers, providing water resources to approximately 1.6 billion people (Bibi et al., 2018). Its unique natural environment and diverse habitats facilitate the interaction and integration of various biota (Chu et al., 2024; Yu et al., 2021). These features help mitigate global warming, protect biodiversity, and provide carbon sequestration, thereby playing a crucial role in ecological security (Liang and Song, 2022). Over the past 50 years, climate change has caused the mean temperature of the TP to rise at twice the global rate (Yang et al., 2019; Zhang et al., 2018a). Plant phenology, which reflects climatic patterns and influences climate through biological and biochemical processes, is a key aspect of this change (Piao et al., 2019; Zhang et al., 2022b). Specifically, the green-up date (D_GU), marking the onset of vegetation growth under favorable climatic conditions (Zhang et al., 2003), is a vital indicator for studying climate change (Shen et al., 2015). Therefore, understanding the changes and mechanisms affecting D_GU is essential for assessing the impact of climate change on the TP and its ecological stability.

Previous studies on climatic factors influencing D_GU have primarily focused on temperature and precipitation, as these are the two most significant determinants of D_GU (Zhang et al., 2022a). Temperature regulates thermal conditions, while precipitation provides the necessary moisture for vegetation growth (Shen et al., 2022). However, the role of snow in influencing D_GU should not be ignored (Wang et al., 2018a). Snow affects vegetation mainly by influencing soil temperature, soil moisture, soil nutrients, and photosynthetically active radiation. Specifically, snow cover can mitigate the exchange of soil heat and moisture with the atmosphere (Zhu et al., 2019a), leading to smaller annual variations in soil temperature compared to air temperature (Zhang et al., 2018b), thereby providing a more favorable overwintering environment for vegetation (Zhu et al., 2019a). Furthermore, snowmelt increases soil moisture, promoting vegetation growth (Peng et al., 2010; Potter, 2020). Changes in soil temperature and moisture due to snow cover can also influence microbial activity, which in turn affects the availability of nutrients for vegetation (Ren et al., 2020; Wang et al., 2015). Additionally, the high albedo of snow reflects much of the solar radiation essential for photosynthesis (Rixen et al., 2022; Yang et al., 2022). Consequently, changes in snow cover invariably influence transitions in D_GU.

The TP has abundant snow cover. In high-altitude areas (around 17.69 %), snow cover typically begins in October and ends in April, whereas in low-elevation areas (around 56.69 %), the snow cover duration is usually less than 20 d, concentrated in December and January (Xu et al., 2024). Snow phenology serves as a crucial indicator of changes in snow cover, and several studies have analyzed the impact of snow phenology on D_GU in the TP. The snow cover end date typically exhibits a significant positive correlation with D_GU, with each 1 d advancement leading to a 0.56 d earlier D_GU (Potter, 2020; Wu et al., 2023). In contrast, the effect of snow cover duration on D_GU is more complex and region-dependent. For instance, a longer snow cover duration leads to a delayed D_GU in the western TP, while it advances D_GU in the eastern TP (Huang et al., 2019; Xiong et al., 2019). Notably, D_GU is most sensitive to the start of snowmelt (D_SOM) among various snow cover phenology metrics on the TP (Xu et al., 2022a). Wang et al. (2015) found that 39.9 % of meadows and 36.7 % of steppes on the TP showed a significant correlation between D_GU and D_SOM. Additionally, Xu et al. (2022a) further highlighted the role of environmental factors in modulating this relationship, noting that precipitation-enhanced humidity strengthened the positive correlation between D_SOM and D_GU, whereas in high-temperature areas, a strong negative correlation prevailed. Besides temperature and precipitation, snowmelt significantly altered environmental conditions conducive to vegetation growth, suggesting that the magnitude of snowmelt between D_SOM and D_GU may also influence ΔD – the time difference between D_GU and D_SOM across the TP (An et al., 2022; Rixen et al., 2022).

Although previous studies have confirmed the role of snowmelt in vegetation green-up dynamics, the response time of vegetation to snowmelt remains poorly understood. This study aims to investigate how long it takes for D_GU to respond to D_SOM and how this response varies regionally. To address this, we calculated and analyzed the time difference between D_GU and D_SOM across the TP, denoted as ΔD. Accurate extraction of D_SOM and D_GU is essential for this analysis. D_SOM detection presents challenges due to the limitations of different remote sensing techniques: while optical remote sensing effectively detects the presence or absence of snow, it struggles to identify the melting state. In contrast, microwave remote sensing provides more reliable snowmelt detection but suffers from lower spatial resolution. To balance these trade-offs, we used a daily snow depth dataset with high spatial resolution to extract D_SOM for the TP from 2001–2018. Similarly, D_GU was derived from an existing high-accuracy dataset for the same period. We then calculated ΔD to assess the delayed response of vegetation green-up to snowmelt. Furthermore, we employed exploratory spatial data analysis to examine the spatiotemporal variation in ΔD, hypothesizing that regional heterogeneity in ΔD is influenced by multiple environmental factors.

2 Materials and methods

2.1 Study area

The TP, located between 25.99–39.82° N and 73.46–104.67° E, has an average elevation exceeding 4000 m (Fig. 1). The TP is characterized by distinct climatic patterns, including intense solar radiation and significant diurnal temperature variations. Winter temperatures range from −15 to −2 °C, while summer temperatures average between 8 and 18 °C. Annual precipitation is approximately 400 mm, with the region transitioning from humid in the southeast to sub-humid, semiarid, and arid conditions in the northwest (Diao et al., 2021).

Figure 1 Map of the study area with vegetation types, distribution of mountains, and geographical zones (base map from ESRI).

Due to substantial topographic uplift and its extensive spatial extent, the TP exhibits a diverse range of climate types and ecosystems. Predominant vegetation types include alpine meadows and alpine steppes, which are widespread across the central regions. The southeastern areas are primarily covered by forests and shrubs, while the western and northern regions are largely barren or desertified due to the terrain and climatic conditions (Zhao et al., 2011).

Snow cover on the TP shows clear spatial and seasonal variability. Except for the Himalayan Mountains, the Pamir Plateau, and the eastern sector of the Tanggula Mountains, which have perennial snow cover, other areas experience seasonal snow cover. Snow accumulation generally begins in September, persists from October to late November, and peaks from December to February. The mean snow cover fraction across the TP is 36.6 % (Zou et al., 2022).

2.2 Data sources

2.2.1 Green-up date dataset

Remote-sensing-based vegetation indices (VIs) are commonly used to identify the D_GU. However, different VIs and extraction methods can introduce uncertainties into the results (Shen et al., 2014). Xu et al. (2022b) compared six different vegetation indices and four extraction methods, finding that the normalized difference greenness index (NDGI) combined with the maximum curvature change rate method (CCRmax) yielded the highest identification accuracy (r=0.62, RMSE=11 d, p<0.01) on the TP. Since the authors of Xu et al. (2022b) are also the collaborators of this study, they directly provided these data. The dataset, spanning 2001–2018, had a spatial resolution of 500 m and was measured in day of year (DOY). To ensure consistent data resolution, we applied bilinear interpolation to reproject the dataset to a resolution of 0.05°.

2.2.2 Snow depth and snowmelt products

Daily snow depth data from 2001–2018 were used to identify the D_SOM. This dataset is accessible through the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/, last access: 20 January 2025) and has a spatial resolution of 0.05°. Produced by Yan et al. (2022), it is based on a long-term series of daily snow depth data in China (1979–2023) (http://poles.tpdc.ac.cn/en/, last access: 20 January 2025) and snow cover probability data derived from the Moderate Resolution Imaging Spectroradiometer (MODIS). Yan et al. (2022) employed a spatial–temporal downscaling method to improve spatial resolution from 0.25° to 0.05°, achieving greater accuracy compared to the original dataset (RMSE=0.61 cm).

Snowmelt data for the period between 2001 and 2018 were sourced from ERA5-Land, provided by the Climate Data Store (https://cds.climate.copernicus.eu/datasets, last access: 20 January 2025). This dataset integrates physical models with global reanalysis observations and has a spatial resolution of 0.1° and a temporal resolution of 1 h. The snowmelt variable quantifies the total amount of water (m) produced by snow melting within snow-covered areas. The data were reprojected to a 0.05° resolution using bilinear interpolation, and the time resolution was adjusted to daily values through summation.

2.2.3 Meteorological dataset

Daily average air temperature and precipitation data from 2001–2018 were used to investigate the impact of meteorological factors on ΔD. These data were obtained from the China Meteorological Forcing Dataset (CMFD), provided by the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/, last access: 20 January 2025). The CMFD integrates remote sensing data with field observations, and due to its use of numerous actual observation sites, it offers higher accuracy compared to the Global Land Data Assimilation System (GLDAS) and ERA5-Land (He et al., 2020; Li et al., 2022). The CMFD dataset has a spatial resolution of 0.1° and a temporal resolution of 3 h, with units of kelvin (K) for temperature and millimeters per hour (mm h⁻¹) for precipitation. To match our analysis requirements, we reprojected the dataset to a 0.05° resolution using bilinear interpolation and adjusted the temporal resolution to daily values by summing (for precipitation) and averaging (for air temperature).

2.2.4 Land cover type

The land cover type map of the TP (2010), which provided vegetation type information for this study, was sourced from the Science Data Bank (https://www.scidb.cn/en, last access: 20 January 2025). This map was classified using the support vector machine method based on MODIS images, achieving a classification accuracy of 93 % with a spatial resolution of 500 m. The map delineates 10 land cover types on the TP: alpine desert, alpine steppe, alpine meadow, bare land/desert, arable land, shrub, coniferous forest, broad-leaved forest, permanent snow/glacier, and lake. Considering that some land covers are non-seasonal or non-vegetation, this study focuses exclusively on alpine steppe, alpine meadow, shrub, coniferous forest, and broad-leaved forest. Prior to use, the dataset was reprojected to a resolution of 0.05°.

2.3 Method

2.3.1 Calculation of ΔD

The time differences between D_GU and D_SOM (ΔD) can be calculated using the following equation:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }D=D_{\text{GU}}-D_{\text{SOM}},\end{array}$$

where D_GU and D_SOM are the green-up date and start of snowmelt, respectively, in DOY.

The D_GU dataset was provided by Xu et al. (2022b). They used the MODIS surface reflectance product MOD09A1 to calculate the NDGI, which integrates red (band 1), near-infrared (band 2), and green (band 4) reflectances. The NDGI time series was then fitted with a four-parameter logistic function, and the curvature change rate (CCR) of the logistic function was calculated. D_GU was defined as the date on which the CCR reached its first local maximum.

D_SOM was identified from snow depth measurements. A snowfall event typically involves two stages: snow accumulation and snowmelt (Fontrodona-Bach et al., 2023), which are characterized by changes in snow depth. The maximum snow depth marks the transition between these two stages (Zheng et al., 2022) (see the orange point in Fig. S1 in the Supplement). To determine D_SOM, the snow depth time series was first smoothed using Savitzky–Golay filtering with a 5 d window to minimize the influence of outliers (as indicated by the orange line in Fig. S1). The date of the maximum snow depth in spring was then identified as D_SOM. Given that multiple snow accumulation and melting events may occur throughout the snow year on the TP (Lei et al., 2023), two criteria were established. First, the snow depth time series was analyzed starting from DOY 61 (1 March), as snowmelt before spring has a minimal impact on vegetation dormancy and snowmelt on the TP typically begins in March (Dong et al., 2024). Second, the snowfall event with the longest duration of snow cover was considered the most significant. Thus, the number of consecutive snow cover days was compared to determine the turning point. According to Stanislaw et al. (2023), a snow cover day is defined as a day when the snow depth exceeds 1 cm. As illustrated in Fig. S1, Snow depth was above 1 cm from DOY 61 to DOY 107, making the duration of the first snowfall event 47 d. In contrast, the durations of the second and third events were 3 and 6 d, respectively. The first snow accumulation and melting process was the longest; thus, the turning point on DOY 96 was designated as D_SOM. Additionally, to ensure that the snow cover is not transient and can influence vegetation, we have introduced a third criterion: the snow cover duration in winter must exceed 10 d (Zhao et al., 2022).

2.3.2 Exploratory spatial data analysis

Classical statistical analysis models assume sample independence, which limits their ability to reveal correlations between the geographical locations of spatial data. Exploratory spatial data analysis (ESDA) offers an enhancement over traditional methods by examining spatial discrepancies and autocorrelation in spatial datasets. Two commonly used indicators in ESDA are the global and local Moran's I. In this study, these indices are applied to analyze the spatial distribution of ΔD. The global Moran's I (I) measures the overall similarity of ΔD between a region and its neighboring regions, as detailed in Eq. (2). The local Moran's I (I_i) measures the degree of clustering or spatial autocorrelation of ΔD within a specific region i, relative to all other regions, as calculated using Eq. (3).

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }I={\displaystyle \frac{n{\sum }_{i=\mathrm{1}}^n{\sum }_{j=\mathrm{1}}^n{w}_{ij}(x_i-\stackrel{\mathrm{‾}}{x})(x_j-\stackrel{\mathrm{‾}}{x})}{{\sum }_{i=\mathrm{1}}^n{\sum }_{j=\mathrm{1}}^n{w}_{ij}{\sum }_{i=\mathrm{1}}^n(x_i-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}}},\text{(3)}& {\displaystyle }& {\displaystyle }{I}_i={\displaystyle \frac{n(x_i-\stackrel{\mathrm{‾}}{x}){\sum }_{j=\mathrm{1}}^n{w}_{ij}(x_j-\stackrel{\mathrm{‾}}{x})}{{\sum }_{i=\mathrm{1}}^n(x_i-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}}},\end{array}$$

where x_i and x_j represent ΔD for region i and j, $\stackrel{\mathrm{‾}}{x}$ is the mean value of ΔD, w_ij is the spatial weight between region i and region j determined by inverse Euclidean distance, and n is the total number of regions.

Moran's I ranges from −1 to 1. A positive value indicates significant spatial autocorrelation, suggesting a tendency for spatial clustering. Conversely, a negative value signifies significant spatial autocorrelation, indicating a tendency for spatial dispersion. If the value is close to zero, the data are considered to be randomly distributed. The significance of spatial autocorrelation can be assessed using the p value of the standardized Z statistic (Eq. 4).

$$\begin{array}{}\text{(4)}& Z={\displaystyle \frac{I-E\left(I\right)}{\sqrt{\text{Var}\left(I\right)}}},\end{array}$$

where I denotes Moran's index, E(I) represents the expectation value of I, and Var(I) is the variance of I. A significance level of p<0.01 is used to determine statistical significance. The calculations were performed using the “Spatial Autocorrelation” and “Cluster and Outlier Analysis” tools in ArcGIS 10.8.

2.3.3 Mann–Kendall test

The Mann–Kendall test was used to investigate the temporal trend of ΔD from 2001–2018 for each pixel. Unlike traditional regression or trend-fitting methods, this test does not require the samples to adhere to a specific distribution. Additionally, the Mann–Kendall test is robust to outliers, thereby minimizing disturbances (Howell et al., 2012).

First, the statistic S must be calculated, as shown in Eq. (5). This statistic represents the sum of the indicator functions that compare the differences between years: if the latter year has a higher value than the former, the function returns 1; if it has a lower value, it returns −1; otherwise, it returns 0 (Semmens and Ramage, 2013). The standardized statistic (Z_c) is then calculated using Eq. (6).

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }S=\sum _{i=\mathrm{1}}^{n-\mathrm{1}}\sum _{j=i+\mathrm{1}}^n\text{sgn}(x_j-x_i)\text{(6)}& {\displaystyle }& {\displaystyle }{Z}_{\mathrm{c}}=\left\{\begin{array}{ll}\frac{S+\mathrm{1}}{\sqrt{\text{Var}\left(S\right)}},& S<\mathrm{0}\\ \mathrm{0},& S=\mathrm{0}\\ \frac{S-\mathrm{1}}{\sqrt{\text{Var}\left(S\right)}},& S>\mathrm{0}\end{array}\right.\end{array}$$

Here n is the length of the time series $x_{\mathrm{1}},\mathrm{\dots },x_i,x_j,\mathrm{\dots }{x}_n$, Var(S) is the variance of the statistic, and Z_c is used for trend testing. If $\left|Z_{\mathrm{c}}\right|>Z_{\mathrm{1}-\mathit{\alpha }/\mathrm{2}}$, a significant temporal trend is present. Here, $Z_{\mathrm{1}-\mathit{\alpha }/\mathrm{2}}$ represents the standard normal variance, and α signifies the significance level. In this study, the significance levels were set at three thresholds: α=0.01, α=0.05, and α=0.1.

If the time series exhibits a significant temporal trend, the trend in ΔD can be assessed using the parameter β (Eq. 7).

$$\begin{array}{}\text{(7)}& \mathit{\beta }=\text{median}\left({\displaystyle \frac{x_i-x_j}{i-j}}\right),\end{array}$$

where x_i and x_j are ΔD of year i and j, respectively. The median function is denoted by “median”, and β indicates the degree of trend: an upward trend if β>0 and a decline if β<0.

2.3.4 The influence of meteorological and snow factors on ΔD

We selected the daily mean temperature of spring (T_spring), the total spring rainfall (P_spring), and the daily mean snowmelt from D_SOM to D_GU (S_StoG) as the influencing factors of ΔD (Shen et al., 2022). The partial correlation coefficient between each variable and ΔD was calculated to quantify their relationship. To isolate the direct effects of T_spring, P_spring, and S_StoG on ΔD while minimizing indirect influences, the second-order partial correlation coefficient was employed (Eqs. 8 and 9).

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{r}_{ij\cdot n}={\displaystyle \frac{r_{ij}-r_{in}\times r_{jn}}{\sqrt{\left(\mathrm{1}-r_{\text{in}}^{\mathrm{2}}\right)\left(\mathrm{1}-r_{\text{jn}}^{\mathrm{2}}\right)}}},\text{(9)}& {\displaystyle }& {\displaystyle }{r}_{ij\cdot mn}={\displaystyle \frac{r_{ij\cdot n}-r_{im\cdot n}\times r_{jm\cdot n}}{\sqrt{\left(\mathrm{1}-r_{im\cdot n}^{\mathrm{2}}\right)\left(\mathrm{1}-r_{jm\cdot n}^{\mathrm{2}}\right)}}},\end{array}$$

where r_ij represents the simple correlation coefficient between variables i and j, calculated using the corr function in Python 3.8. r_ij⋅n is the first-order partial correlation coefficient between i and j, which accounts for their correlation after removing the linear effects of control variable n. Similarly, r_ij⋅mn is the second-order partial correlation coefficient between i and j, controlling for both m and n. Accordingly, we calculated $r_{\mathrm{\Delta }D{T}_{\text{spring}}\cdot P_{\text{spring}}{S}_{\text{StoG}}}$, $r_{\mathrm{\Delta }D{P}_{\text{spring}}\cdot T_{\text{spring}}{S}_{\text{StoG}}}$, and $r_{\mathrm{\Delta }D{S}_{\text{StoG}}\cdot T_{\text{spring}}{P}_{\text{spring}}}$.

Next, a multiple linear regression model was established for each pixel (Eq. 10). The partial correlation coefficient between each variable and ΔD was calculated to quantify their relationship. Subsequently, a multiple linear regression model was established for each pixel (Eq. 8). A prerequisite for multiple linear regression is passing the collinearity test, which requires the variance inflation factor (VIF) to be less than 3, indicating no collinearity. In this study, the VIF values for T_spring, P_spring, and S_StoG were 1.243, 1.209, and 1.174, respectively, confirming that the collinearity test was satisfied. The OLS and variance_inflation_factor functions were used to construct the linear regression model and perform the collinearity test in Python 3.8.

$$\begin{array}{}\text{(10)}& \mathrm{\Delta }D=a{T}_{\text{spring}}+b{P}_{\text{spring}}+c{S}_{\text{StoG}}+d,\end{array}$$

where a,b, and c are the fitting coefficients used to determine the contribution of T_spring, P_spring, and S_StoG to ΔD, respectively, while d represents a constant term. By comparing the fitting coefficients, the factor with the largest coefficient was deemed dominant.

It is important to note that all variables were standardized by Z scores before being input into the model to eliminate dimensional differences. Given that the model was based on the pixel scale and the sample size was too small, we have expanded the sample pixel to its surrounding eight neighboring pixels. The significance level was set at p<0.05.

Figure 2 Spatial distribution of multiyear-averaged (a) D_SOM and (b) D_GU from 2001–2018 on the Tibetan Plateau.

3 Result

3.1 Spatial variation in D_SOM and D_GU

Figure 2 shows the spatial pattern of the multiyear mean D_SOM and D_GU data from 2001–2018. The Qaidam Basin and southern Qiangtang Plateau, which rarely experience snowfall (Xu et al., 2024), present challenges for snowmelt detection. In the northern TP, where vegetation is sparse, and in the southeastern TP, where seasonal vegetation changes are minimal, D_GU was not observed. Generally, D_SOM occurred between DOY 70–100 (10 March–9 April) for 73.5 % of the TP, with an average DOY of 86.7. The latest D_SOM (>DOY 105) was found in the Kunlun and Nyainqentanglha mountains, while the earliest D_SOM occurred in the Hengduan and Minshan mountains, with an averaged DOY of 76.8. In contrast, D_GU appeared later, ranging from DOY 110–140 (19 April–19 May) for 77.6 % of the TP. Vegetation greens up earlier in the southeast and later in the center and west.

Figure 3 (a) Spatial and frequency distribution histogram and (b) local and global Moran's I values of average ΔD on the Tibetan Plateau over 2001–2018.

3.2 Spatial and temporal distribution characteristics of ΔD

3.2.1 Spatial distribution of ΔD

ΔD was calculated using Eq. (1) for regions where both D_SOM and D_GU were non-null (Figs. 3a and S2 in the Supplement). A positive ΔD indicates that D_SOM preceded D_GU, and vice versa. The average ΔD over the TP from 2001–2018 was 38.5 d, with positive values observed over 95 % of the TP. The proportion of ΔD ranging from 40–60 d was the highest, accounting for 40.4 % on average. In the east TP, ΔD increased from south to north, with larger values observed in the Gangdise and Himalayas mountains. Negative ΔD was mainly found in the northwest TP.

To gain more insight about spatial distribution characteristics of ΔD, the global and local Moran's I were calculated (Figs. 3b, S3, and S4 in the Supplement), The global Moran's I of ΔD over the TP from 2001–2018 was 0.3 and was significant at the 0.01 level ($p=\mathrm{0.00}<\mathrm{0.01}$). This indicates that ΔD demonstrates spatial clustering. The local Moran's I index revealed that about 81.7 % of the regions exhibited significant spatial autocorrelation. The pale blue regions in Fig. 3b represent low-value clustering, while the pink regions indicate high-value clustering, meaning ΔD was generally longer in these areas. Red (blue) areas represent high (low) values surrounded by low (high) values. Low values are predominantly concentrated in the Hengduan Mountains, and high values are primarily clustered around the Bayankara and Himalaya mountains. In the eastern TP, the Bayankara Mountains mark the boundary between high and low values. In the central and northwestern TP, negative ΔD results in the low-value cluster with scattered positive values.

3.2.2 Temporal dynamic of ΔD

Regarding the temporal trend of ΔD, Fig. 4 shows the interannual variation from 2001–2018, along with its significance indicated in the bottom-left corner. Only 5.2 % of the study area, predominantly in the eastern TP, exhibits a significant trend. Of these regions, approximately 36.6 % show a declining trajectory, suggesting that vegetation in these areas responded more rapidly to snowmelt. However, for most regions, no discernible trend was observed over time, likely because the warming hiatus since the 2000s has not led to noticeable advancements in D_SOM and D_GU (Piao et al., 2019; Wang et al., 2019).

Figure 4 Interannual variation trend and significant levels of ΔD on the Tibetan Plateau from 2001–2018.

3.3 Influence of air temperature, precipitation, and snowmelt on ΔD

As observed in Sect. 3.2, ΔD is positive in most regions. Furthermore, snowmelt only significantly impacts vegetation when D_SOM precedes D_GU. Consequently, subsequent analyses will focus exclusively on areas where ΔD is positive.

Figure 5 illustrates the mean value of ΔD under varying spring meteorological conditions. ΔD exhibits a clear stepwise decline from cold to warm regions, decreasing from approximately 48 to 37 d (Fig. 5a). In colder or hotter spring conditions (i.e., T_spring<270 K or T_spring>275 K), ΔD decreased slightly. However, near the freezing point (270–275 K), ΔD shortens by 3 d with each 1 K increase in T_spring. Under various precipitation conditions (Fig. 5b), ΔD shortens by 0.29–1.96 d for every 10 mm increase in P_spring. Figure 5c reveals a strong negative correlation between ΔD and S_StoG when S_StoG exceeded 6 mm d⁻¹. For each 1 mm increase in S_StoG, ΔD decreases by approximately 0.615 d. The dispersion within each snowmelt category remains relatively consistent, with a standard deviation of about 16.8 d.

Figure 5 Variations in ΔD across regions with differing (a) spring mean temperature (T_spring), (b) spring total rainfall (P_spring), and (c) daily snowmelt from D_SOM to D_GU (S_StoG). Points represent the mean ΔD, while error bars denote 1 standard deviation. The slope and R² value reflect the coefficient and precision of the linear regression, respectively, with a significance level of 0.01.

To quantify the influence of environmental factors on ΔD across different regions, a partial correlation coefficient was calculated for each pixel (Fig. 6). Only pixels with valid ΔD over a 6-year period were included in the analysis.

Figure 6 Spatial distribution of the partial correlation between ΔD and (a) spring mean temperature ($R_{\mathrm{\Delta }D\text{ and }{T}_{\text{spring}}}$), (b) spring total rainfall ($R_{\mathrm{\Delta }D\text{ and }{P}_{\text{spring}}}$), and (c) daily snowmelt from D_SOM to D_GU ($R_{\mathrm{\Delta }D\text{ and }{S}_{\text{StoG}}}$) with a significance level of 0.05.

At a significance level of 0.05, ΔD was significantly correlated with temperature in 23.5 % of the samples. Among these significant pixels, 51.3 % exhibited a positive correlation between ΔD and T_spring. These positively correlated pixels were predominantly located in the southern part of the valid data, while negatively correlated pixels were concentrated in the northern regions. Regarding rainfall, 28.8 % of the samples exhibited a significant correlation with ΔD, with 82.9 % of these showing a negative correlation and 17.1 % showing a positive one. Snowmelt demonstrated a significant relationship with ΔD in 35.4 % of samples, with 83.8 % showing a negative correlation, and the average correlation index was −0.27.

A linear regression model was established to determine the dominant factor influencing ΔD at each pixel (Figs. 7 and S5 in the Supplement). Among the significant samples (p<0.05), T_spring, P_spring, and S_StoG accounted for 34.1 %, 20.3 %, and 45.6 % of the dominant factors, respectively. On the southeastern TP, T_spring primarily influences ΔD, whereas on the northwestern TP, S_StoG is more influential. This aligns with the distribution of vegetation types and geographical zones (Fig. 1). Consequently, the proportion of each dominant factor across different regions was assessed. In the arid zone, 67 % of pixels were dominated by S_StoG, compared to 20 % dominated by T_spring. As moisture increased, the dominant factor shifted from S_StoG to T_spring, with 48 % of pixels in humid regions under temperature dominance and only 18 % under S_StoG. Similarly, the proportion of T_spring-dominant pixels increased from 20 % to 44 %, while S_StoG-dominant pixels decreased from 64 % to 27 % as vegetation became taller. Additionally, P_spring-dominant pixels consistently represented the lowest fractions (12 %–36 %) across all regions.

Figure 7 (a) Spatial distribution of dominant factor of ΔD with a significance level of 0.05 and its proportion diagram among (b) different geographical zones and (c) different vegetation types.

4 Discussion

4.1 Quality evaluation of ΔD

ΔD is primarily determined by D_SOM and D_GU. Xu et al. (2022b) discussed the accuracy of D_GU, so this study focused on the identification quality of D_SOM.

D_SOM can currently be derived from optical or microwave remote sensing images. The normalized difference snow index (NDSI) is commonly used for snow identification in optical images. During periods of surface snow coverage, the NDSI tended to remain elevated and steady. After the start of snowmelt, the NDSI began to decrease. Consequently, D_SOM can be identified using the dynamic threshold method on the NDSI time-series curve (Potter, 2020; Zheng et al., 2022). However, NDSI mainly reflects snow presence rather than the snowmelt process, which initially reduces snow depth rather than snow cover extent (Panday et al., 2011).

Microwave bands, on the other hand, are more effective at tracking the internal state of snowpacks due to differences in the dielectric constant between water and snow (Ma et al., 2020). Melting snow leads to abrupt changes in brightness temperature or backscatter coefficient, facilitating D_SOM identification (Smith et al., 2017). Methods developed for this purpose include the band threshold method (Howell et al., 2012), the diurnal amplitude variation algorithm (DAV) (Semmens and Ramage, 2013), and the cross-polarized gradient ratio (Grippa et al., 2005). For instance, Xiong et al. (2017) identified D_SOM as the date when the differential average derivative of the Ku-band reached a maximum, providing the snowmelt onset time for the SOMHMA high-mountain Asia (1979–2018) dataset, available from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/, last access: 20 January 2025) with a spatial resolution of 0.25°.

Although this dataset captures snowmelt evolution, its spatial resolution is relatively low. D_SOM can also be indirectly identified from abrupt points in the snow depth/snow water equivalent time series (Fontrodona-Bach et al., 2023; Zheng et al., 2022). Our study used the snow depth dataset from Sect. 2.2.2, which improved spatial resolution to 0.05° via downscaling while maintaining the benefits of microwave bands, making it more accurate than traditional microwave data with a 0.25° resolution.

Verifying D_SOM is challenging due to the lack of in situ data. Snowmelt patterns, influenced by external conditions and similar to temperature distributions, have been cross-validated with temperature records in previous studies (Drobot and Anderson, 2001; Panday et al., 2011; Zheng et al., 2020). Thus, following Grippa et al. (2005), we randomly selected one sample area for each land cover on the TP and calculated the correlation coefficient between average D_SOM and mean April temperature (Fig. 8). Our results showed stronger consistency with temperature trends compared to SOMHMA. The ΔD calculated by SOMHMA was relatively random and inconclusive (Figs. S6–S8 in the Supplement), potentially due to mixed-pixel issues from its limited resolution. Although our method offers improved accuracy, some errors in D_SOM persist due to original data uncertainty (RMSE=0.61 cm) (Yan et al., 2022).

Figure 8 Relevance between the average D_SOM and mean April temperature in a sample area of each land cover type.

4.2 Response of vegetation to snowmelt

In this study, the time difference between D_GU and D_SOM was defined as ΔD. This metric can have both positive and negative values, each with different implications. A positive ΔD indicates the time it takes for vegetation to respond to snowmelt. On the TP, ΔD was predominantly positive, with an average of 38.5 d. This suggests that vegetation response to snowmelt typically involves a delay of nearly 1 month. In contrast, response times are generally shorter in other regions. ΔD in the high Arctic, Alaska, and Finland averaged 3 weeks (Cooper et al., 2011), 14.9 d (Zheng et al., 2022), and 5–13 d (Manninen et al., 2019), respectively. They were shorter than those on the TP, likely due to differences in snowpack characteristics. The blocking effect of the Himalayas and Karakoram mountains results in transient, shallow, and patchy snow on the TP (Lei et al., 2023). Consequently, the impact of snow on response times is particularly significant on the TP.

The influence of S_StoG on response rate was notably widespread, with more S_StoG shortening the response time in one-fourth of the TP. It can be inferred that reductions in ice volume and increases in the rain-to-snow ratio due to global warming may influence vegetation sensitivity to snowmelt (Tomaszewska et al., 2020). Snowmelt primarily affects soil moisture and nutrients, thus stimulating vegetative response to D_SOM. Previous studies have emphasized the critical role of snowmelt in soil moisture. During early vegetation development, especially in alpine and cold regions, snow serves as an essential water source (Ernakovich et al., 2014; Li et al., 2020; Ma et al., 2024). This is evident in Fig. 7, where regions dominated by S_StoG were more extensive than those dominated by P_spring. Edwards et al. (2007) observed that nitrogen is leached from the snow layer during snowmelt and combines with nitrogen released from mineralization to form a pulse of inorganic nitrogen beneath the snow. This results in a peak in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ available post-snowmelt. As alpine vegetation typically prefers ${\mathrm{NO}}_{\mathrm{3}}^{-}$ as a nitrogen source, snowmelt meets this need (Broadbent et al., 2021). Additionally, the insulating effect of snow on vegetation diminishes after melting (Rixen et al., 2022). For instance, Starr and Oberbauer (2003) measured solar radiation at different snow depths in the northern foothills of the Brooks Mountain Range, Alaska, finding that light levels were reduced by 20 %–40 % under 10 cm of snow and declined by up to 70 % under 20 cm compared to bare ground. Thus, future remote sensing efforts should aim to quantify the snow's shielding effect to better understand vegetation response mechanisms.

The vegetation response to D_SOM is influenced not only by snowmelt but also by meteorological conditions. Temperature and precipitation, as primary factors of D_GU, significantly affect the response rate on the TP, accounting for 23.5 % and 28.8 %, respectively (Fig. 6). Precipitation provides essential water for vegetation, enabling a swift return to growth post-snowmelt. Temperature, on the other hand, impacts the heat required for vegetation growth and accelerates snowmelt, thereby supplying water and nutrients (Liu et al., 2021). Higher temperature can thus enhance D_GU and shorten response time (Fig. 5a). However, elevated temperatures may also advance D_SOM (Mioduszewski et al., 2014; Mioduszewski et al., 2015). Due to the differential effects of temperature on D_SOM and D_GU, both a positive and a negative correlation between T_spring and ΔD were observed at local scales (Fig. 6a). Statistical analysis shows that positively correlated pixels are typically found in southern regions with an average T_spring of 271.19 K, while negatively correlated areas are in northern regions with an average T_spring of 269.84 K. In warmer regions with mean annual temperatures above freezing, spring temperature correlates negatively with D_SOM ($r=-\mathrm{0.46}$) and D_GU ($r=-\mathrm{0.07}$), indicating that temperature primarily influences snowmelt rather than vegetation growth, thus extending response times. In colder regions increased temperatures can reduce cold stress on vegetation, resulting in a larger effect on D_GU ($r=-\mathrm{0.27}$). However, consistent sub-freezing temperatures do not significantly lead to later $D_{\text{SOM}}(r=-\mathrm{0.28}$, which is similar with D_GU). In summary, the relationship between temperature and response time is modulated by the magnitudes of their respective influences at the local scale.

Furthermore, the primary factors controlling response time varied across regions but exhibited a clear regional distribution. Previous studies have indicated that in arid areas, D_GU is more sensitive to water availability, especially in the central and western TP, while D_GU shows a more significant negative correlation with temperature in the eastern TP, which experiences wetter springs (Piao et al., 2019; Shen et al., 2022). Despite differences in dependent variables, this conclusion remains consistent (Fig. 7b). In arid regions, where water is the main limiting factor, ΔD is predominantly controlled by snowmelt. In contrast, in humid regions where heat is the limiting factor, ΔD is mainly influenced by T_spring. When analyzing by vegetation type, T_spring was the dominant factor for higher vegetation (e.g., shrubs and forests), while lower vegetation (e.g., alpine steppes and meadows) is more affected by S_StoG (Fig. 7c). Wind speed is high on the TP due to its special topography and elevation. This limited the snow accumulation on the branches for higher vegetation, and thus some parts of the vegetation are exposed to the outside environment. In contrast, the height of the lower vegetation is sometimes equal to the depth of the snow, which can partially or completely cover the vegetation (Tang et al., 2024). If the temperature rises, it will directly act on higher vegetation while heat will transfer slower to the low vegetation buffered by snow. If the snow melts, the melting water can directly affect the lower vegetation under the snow while the tall vegetation can only be used after absorption through the soil. This finding aligns with the observations of Zheng et al. (2022), who noted that soil temperature increased more rapidly in lower vegetation than in higher vegetation following snowmelt in Alaska, suggesting that snowmelt has a greater impact on lower vegetation. Our study, which compares these factors at the pixel scale, further substantiates this view and clarifies the spatial variations across the TP.

Specially, if ΔD is negative (i.e., D_SOM occurs later than D_GU), this situation is termed “false spring” (Chamberlain et al., 2019, 2021; Peterson and Abatzoglou, 2014). During false springs, once plant growth begins, frost tolerance significantly decreases. Continued exposure to freezing conditions can damage vegetation tissues and reduce productivity, leading to both ecological and economic impacts (Chamberlain et al., 2019). Zhu et al. (2019b) used the Spring Indices (SI-x) model to simulate the probability of false spring in China from 1950–2005 and projected it until 2100. Their findings revealed that the central TP exhibited the highest probability of false spring in China. Rising temperatures due to global warming are expected to further advance the D_GU, thereby increasing the risk of false spring. Thus, false spring warrants further investigation in the TP, even though it was not the primary focus of this study.

Although previous studies have established a strong positive correlation between D_SOM and D_GU (An et al., 2022; Wang et al., 2015, 2018b; Xu et al., 2022a), the specific time differences between them in the TP remain unclear. Our study advances this understanding by examining the heterogeneity and mechanisms of vegetation response to snowmelt and reinforces the notion that snowmelt primarily affects arid regions and areas with low vegetation cover.

5 Conclusion

This study investigates the dynamic response of vegetation to snowmelt on the Tibetan Plateau from 2001–2018. Our results reveal that the effect of snowmelt on vegetation is not immediate, with a mean response lag of 38.5 d from D_SOM to D_GU. Notably, the false spring was observed in the northwestern TP, which warrants further exploration. As precipitation and snowmelt increase, the response time shortens. More complex than these factors, temperature exerts a greater influence on D_GU than D_SOM in colder regions, thus shortening the response time. Conversely, in warmer areas, the increased temperature has a stronger impact on D_SOM, which lengthens the response time. Furthermore, vegetation in arid regions is more dependent on water than heat, and low-vegetation areas rely more on sub-snow habitats than external climatic factors. These findings provide valuable insights into how vegetation responds to snowmelt in the context of climate change, deepening our understanding of the relationship between snowmelt onset and green-up dates. This knowledge is essential for predicting vegetation phenology and managing ecosystem services under changing climate conditions. Future research should focus on the impacts of snow cover and false spring.