Temperature-dominated spatiotemporal variability in snow phenology on the Tibetan Plateau from 2002 to 2022

A detailed understanding of snow cover and its possible feedback on climate change on the Tibetan Plateau (TP) is of great importance. However, spatiotemporal variability in snow phenology (SP) and its influencing factors on the TP remain unclear. Based on the daily gap-free snow cover product (HMRFS-TP) with 500 m resolution, this study investigated the spatiotemporal variability in snow cover days (SCDs), snow onset date (SOD), and snow end date (SED) on the TP from 2002 to 2022. A structural equation model was used to quantify the direct and indirect effects of meteorological factors, geographical location, topography, and vegetation greenness on SP. The results indicate that the spatial distribution of SP on the TP was extremely uneven and exhibited temporal heterogeneity. SP showed vertical zonality influenced by elevation (longer SCD, earlier SOD, and later SED at higher elevations). A total of 4.62 % of the TP area had a significant decrease in SCDs, at a rate of −1.74 d yr⁻¹. The SOD of 2.34 % of the TP area showed a significant delayed trend, at a rate of 2.90 d yr⁻¹, while the SED of 1.52 % of the TP area had a significant advanced trend, at a rate of at −2.49 d yr⁻¹. We also found a strong elevation dependence for the trend in SCDs ($R=-\mathrm{0.73}$). Air temperature, precipitation, wind speed, and shortwave radiation can directly affect SP as well as indirectly affect it by influencing the growth of vegetation, whereas the direct effect was much greater than the indirect effect. Geographical location (latitude and longitude) and topographic conditions (elevation and slope) indirectly affected SP by modulating meteorological conditions and the growth of vegetation. Vegetation primarily influences SP by intercepting the snow and regulating the balance of the solar radiation budget. Regarding the total effect, air temperature was found to be the dominant factor. This study contributes to the understanding of snow variation in response to global warming over the past 2 decades by providing a basis for predicting future environmental and climate changes and their impacts on the TP.

1 Introduction

The rapid accumulation and melting of snow make it one of the most active natural materials on the Earth's surface (Gutzler and Rosen, 1992; Ma et al., 2023). The Tibetan Plateau (TP) is one of the most sensitive regions to climate change because of its unstable alpine ecosystems and fragile natural environment (Huang et al., 2017). In the past 50 years, temperature has risen at a rate of about 0.37 °C per decade on the TP (You et al., 2021). Rising temperatures inevitably affect snow accumulation and melting processes, exemplified by variations in snow phenology (SP), including snow cover days (SCDs), snow onset date (SOD), and snow end date (SED) (Chen et al., 2015; Ma et al., 2020). Variations in SP can in turn affect terrestrial ecosystems and feedback on regional climate (Cherkauer and Sinha, 2008). Early snowmelt also causes substantial variability in the onset date and amount of snowmelt runoff, increasing the incidence of disasters (e.g., floods) (Fyfe et al., 2017). A detailed understanding of the variability in SP and its possible feedback on climate change is of great significance to the hydrological cycle (Kraaijenbrink et al., 2021), ecological balance (Keyser et al., 2022), and societal security (Wang et al., 2017) on the TP.

There are typical interannual variations in the SP on the TP. Ground-measured measurements and remote sensing data have been extensively used in recent decades to better understand these changes (Ma et al., 2022; Notarnicola, 2020). Remote sensing is beneficial in high-elevation areas where few ground-measured measurements are obtainable (Huang et al., 2022a). Passive microwave satellite (e.g., SMMR, SSM/I, and AMSE-R) and multispectral satellite (e.g., Landsat, Sentinel-2, and MODIS) data are widely used to retrieve snow information and reveal its variability (Chen et al., 2018). Although passive microwave data are available over an extended period, the spatial resolution (∼ 10–25 km) is relatively coarse (Huang et al., 2017). Landsat and Sentinel-2 have a high spatial resolution (∼ 10–30 m) but a relatively coarse temporal resolution (10 or 16 d), inadequate for monitoring temporal variability in SP (Huang et al., 2022a). Currently, the most widely available snow data used for the TP are the daily MODIS snow product (Hall et al., 1995), which has revealed that the SP on the TP has undergone changes in the last 2 decades (Ma et al., 2023; Tang et al., 2022; Wang et al., 2017). In addition, owing to the complex topography of the TP, elevation may have a particular impact on SP. However, it is still unclear whether trends in SP are elevation-dependent (Ma et al., 2023; Wang et al., 2021). Moreover, since the original daily MODIS images have data gaps due to cloud cover and relatively low estimation accuracy due to the complex topography of the TP (Huang et al., 2022b), most previous studies used relatively simple methods to fill data gaps and did not consider topographic effects. Using snow products without topographic corrections may result in inaccurate extraction of SP, further affecting the analysis result of interannual variations. Huang et al. (2022b) generated a daily gap-free snow product on the TP by employing a hidden Markov random field (HMRF) framework to the MODIS product (HMRFS-TP). The data gaps were filled by optimally integrating spatiotemporal, spectral, and environmental information. The accuracy of HMRFS-TP was significantly improved over the complex topography, providing spatiotemporally continuous information on snow distribution over an extended period with high accuracy.

Regarding the influencing factors on SP of the TP, previous studies have highlighted a complex and heterogeneous situation in which SP is driven by meteorological factors, topographic conditions, geographical location, and vegetation greenness (Pulliainen et al., 2020; Qi et al., 2021; You et al., 2020). The impact of air temperature and precipitation on snow cover has been extensively studied (Li et al., 2022; Tarca et al., 2022). However, SP is also closely associated with other meteorological factors, such as wind and humidity. But their quantitative influence on SP on the TP has not been discussed further (Ma et al., 2022; Xie et al., 2017). Topographic characteristics (Guo et al., 2022) and vegetation greenness (Qi et al., 2021) are also considered responsible for SP. Elevation is an essential factor that influences SP because snow accumulates more at higher elevations. Steeper slopes may receive less snow than flatter ones due to gravitational pull (Li et al., 2011). Vegetation greenness changes the snow redistribution process by intercepting snow and affects the spatial pattern and melt rate of surface snow by regulating the solar radiation budget balance (Barrere et al., 2018).

Although much research has examined the variability in the SP and its influencing factors on the TP, uncertainties and limitations remain. Previous studies have mainly focused on the response of SP to temperature and precipitation, but few studies have quantitatively explored other associated factors influencing SP. Additionally, many analyses have primarily focused on examining the relationship between SP and individual influencing factors, without taking into account the interrelationships among these factors. Associated factors may have both direct and indirect impacts on SP. For example, temperature can affect SP directly, as well as indirectly by influencing the growth of vegetation (Qi et al., 2021). The specific process by which each factor influences SP (i.e., directly or indirectly) has not been elucidated. The structural equation model (SEM) is a suitable method for studying such a coupling relationship, which can quantitatively explain the direct and indirect relationships between factors (Grace et al., 2010; Shen et al., 2022; Zhang et al., 2022). Here we adopt this approach to explore the mediating effects between SP and associated factors.

In this study, we conducted a comprehensive investigation on the spatiotemporal variability in SP on the TP over the past 2 decades (2002–2022), based on the daily gap-free HMRFS-TP dataset, and further quantified the influence of associated factors on SP. The main goals are to (1) calculate SP parameters (SCD, SOD, and SED) and analyze their spatiotemporal variability, (2) examine whether the trends in SP are elevation-dependent, and (3) quantify the direct and indirect effects of associated factors on SP.

2 Study area and data

2.1 Study area

The TP (26°00^′12^′′ N–39°46^′50^′′ N, 73°18^′52^′′ E–104°46^′59^′′ E; Fig. 1) is the largest snow cover area in the middle latitudes of the Northern Hemisphere, with 10⁵ km² of glaciers and an annual snowfall of 41.9 × 10⁹ m³ (Yao et al., 2012). The southeastern TP is warm and moist, while the northwestern TP is cold and dry. The average annual temperature ranges from −6 to 20 °C and the average annual precipitation ranges from 150 to 800 mm (Wang et al., 2023). Temperature in its permafrost region is rising at about 3 times the speed of global warming (Wang et al., 2022). As a topographic barrier, the TP forces air masses up and produces cooling effects (Wu et al., 2015) and further influences the thermodynamic properties of atmospheric circulation, strength, and duration of Asian monsoon systems, global climate, and energy budget (Fan et al., 2019). It is the source of several prominent Asian rivers (e.g., the Yellow River, Yangtze, Indus, and Mekong), which are primarily fed and regulated by meltwater from glaciers and snow (Chen et al., 2022).

Figure 1 Topography, lakes, and distribution of selected meteorological stations for validation on the Tibetan Plateau (TP) (base map from ESRI).

2.2 Datasets

2.2.1 Snow cover product

The daily gap-free HMRFS-TP at 500 m resolution from 2002 to 2022 was used to extract SP (Huang et al., 2022b). The HMRF technique produces a dataset (Huang et al., 2018) that optimally couples spatiotemporal, spectral, and environmental information to fill data gaps in original MODIS images. The accuracy of the HMRFS-TP is 91.36 % and 98.29 % based on snow maps derived from Landsat images and in situ observations, respectively. The estimation accuracy is notably improved in snow transitional periods and complex topography with higher elevations and sunny conditions. To determine the SP (including SCD, SOD, and SED), we specify the snow season from 1 September of the previous year to 31 August of the current year (Tang et al., 2022). The snow accumulation season is defined as from 1 September of the previous year to 28 (or 29) February of the current year, and the snowmelt season is from 1 March to 31 August of the current year (Ma et al., 2023). The SP for 20 snow seasons from 1 September 2002 to 31 August 2022 was extracted based on the long-term HMRFS-TP dataset.

2.2.2 Meteorological data

We used a daily meteorological dataset at $\mathrm{1}/\mathrm{30}$° resolution from 2002 to 2022 to investigate the impact of meteorological factors on SP, including air temperature, precipitation, specific humidity, wind speed, and downward shortwave radiation. This dataset was generated by integrating in situ observations, remote sensing, and a reanalysis dataset (He et al., 2020; Yang et al., 2023). All daily meteorological data were resampled to 500 m to maintain consistency with the MODIS snow product. We calculated average air temperature, average humidity, average wind speed, total precipitation, and total shortwave radiation for each snow season, snow accumulation season, and snowmelt season to examine their relationship with SCD, SOD, and SED, respectively (Chen et al., 2018; Ma et al., 2023).

2.2.3 Digital elevation model (DEM)

The 30 m DEM data, available from the USGS Earth Explorer (https://earthexplorer.usgs.gov, last access: 21 March 2023), were resampled to 500 m resolution. Topographical parameters including elevation, slope, and aspect were calculated from it.

2.2.4 Land surface reflectance product

Compared with other vegetation indexes, the normalized difference greenness index (NDGI) (Yang et al., 2019) has been proven to more accurately represent the vegetation greenness growth status in snow-covered areas such as the TP (Xu et al., 2022a, b). Here, we calculated the NDGI based on MODIS Terra surface reflectance MOD09A1 with 500 m resolution and an 8 d repeat cycle from 2002 to 2022 (Vermote, 2021). The bands 4, 1, and 2 of MOD09A1 were used to calculate the NDGI. The maximum complex method was implemented to filter noise and fill data gaps owing to cloud and atmosphere in NDGI series (Xu et al., 2022b). The average NDGI was calculated for each snow season, snow accumulation season, and snowmelt season.

2.2.5 Auxiliary data

Snow depth recorded from meteorological stations was used as the reference value (ground-observed SP) to validate the accuracy of the satellite-derived SP parameters. Due to instrument machine failure or human-made errors, not all meteorological stations (approximately 137) on the TP recorded snow depth values every day. This resulted in many data gaps during the snow season, which were not sufficiently continuous to extract ground-observed SP. Therefore, 24 meteorological stations with daily records were selected for validation (Fig. 1).

3 Methods

Figure 2 illustrates the workflow of this study. The SP parameters (SCD, SOD, and SED) were first extracted from the daily HMRFS-TP dataset, and their spatiotemporal variability were analyzed. We then adopted SEM to quantify the direct and indirect effects of associated factors on SP. The impact of relatively important factors on SP was discussed in further depth.

Figure 2 Flowchart to reveal the spatiotemporal variability in snow phenology (SP) and the direct and indirect effects of associated factors on the TP from 2002 to 2022. SEM refers to the structural equation model.

3.1 Snow phenology extraction and accuracy assessment

Three parameters, SCD, SOD, and SED, were derived to describe SP in each snow season and were calculated from the daily HMRFS-TP dataset at the pixel level. SCD is the sum of days when a pixel is covered by snow across the snow season. SOD is the first date when a pixel is covered with snow and stays at least 5 d in a snow season. SED is the last date on which a pixel was identified as snow for 5 consecutive days (Tang et al., 2022). Using a threshold value of 5 d can reduce the influence of frequent short-term variability (e.g., snowmelt and accumulation in early spring and late autumn), which has been a widely used threshold on the extraction of SOD and SED (Guo et al., 2022; Wang et al., 2017; Xu et al., 2022b). Calculations of SP were excluded from regions containing lakes.

The accuracy of the extracted SP was validated using ground-observed SP from snow depth measured by meteorological stations. The ground-observed SP was calculated as follows. First, all recorded snow depth values were reclassified as snow or no snow based on a 3 cm threshold (Huang et al., 2022a, b). For each snow season, meteorological stations that recorded more than 200 d of snow depth were selected (Hao et al., 2022; Zhao et al., 2022). Days without snow depth records at the selected stations were supplemented with remote sensing images (Landsat series, Sentinel-2, etc.) of these days. Stations with fewer than 20 snow-covered days and fewer than 5 consecutive snow-covered days during the snow season were excluded. After applying these criteria, a total of 56 ground-observed SP instances from 24 stations were used for accuracy validation.

3.2 Trend analysis

To explore the interannual variation trend and trend significance in SP from 2002 to 2022, the Theil–Sen nonparametric regression (Sen, 1968) and Mann–Kendall (M–K) tests (Hirsch et al., 1982) were applied to each pixel. The Theil–Sen method has the advantage of dealing with non-normally distributed data and is robust against outliers compared with traditional linear regression (Theil, 1992).

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

where β is the trend slope; a positive value denotes the trend of SP as a delay/extension, and a negative value assumes SP is advanced/reduction. n=20, and x_i is the ith value in 20 years, the same as x_j and jth.

The M–K test is a nonparametric approach for monotonic trend that has been used for trend detection of hydrological and meteorological time series (Qi et al., 2021). The Z value assumes the temporal trend is statistically significant:

$$\begin{array}{}\text{(2)}& {\displaystyle }Z=\left\{\begin{array}{l}\frac{S-\mathrm{1}}{\sqrt{\mathrm{Var}\left(S\right)}},\text{if}S>\mathrm{0}\\ \mathrm{0},\text{if}S=\mathrm{0}\\ \frac{S+\mathrm{1}}{\sqrt{\mathrm{Var}\left(S\right)}},\text{if}S<\mathrm{0}\end{array},\right.\mathrm{Var}\left(S\right)={\displaystyle \frac{n(n-\mathrm{1})(\mathrm{2}n+\mathrm{5})}{\mathrm{18}}},\text{(3)}& {\displaystyle }S={\sum }_{i=\mathrm{1}}^{n-\mathrm{1}}{\sum }_{j=i+\mathrm{1}}^n\left\{\begin{array}{l}\mathrm{1},\text{if}{x}_j>x_i\\ \mathrm{0},\text{if}{x}_j=x_i,\mathrm{1}<i<j<n\\ -\mathrm{1},\text{if}{x}_j<x_i\end{array}\right..\end{array}$$

When |Z|>1.28, 1.64, and 2.32, the tests were significant at levels of 0.1, 0.05, and 0.01, respectively.

Since snow cover may not be present each year, to avoid unnatural results caused by a small amount of data, only the pixels with at least 6 recorded years for the SP were used for the interannual trend evaluation (Xu et al., 2022b).

3.3 Structural equation model

We used SEM to quantify the direct and indirect effects of various factors on SP. It is a multivariate collection of methods that can simulate the interaction between various factors at the same time, supplying a framework for extrapolating cause–effect relationships and revealing direct and indirect relations between independent and dependent variables (Grace et al., 2010). SEM includes a covariance-based structural equation model (CB-SEM) and the partial least-square structural equation model (PLS-SEM) (Venturini and Mehmetoglu, 2019). Contrasted with CB-SEM, PLS-SEM focuses on mining sample information and can reflect the nature and structural characteristics of objects as much as possible, making it more suitable for exploring newly constructed structural models (Hair et al., 2011; Ringle et al., 2012). More importantly, PLS-SEM is a nonparametric model that does not require a normal distribution of samples and has considerable potential for remote sensing applications (Lopatin et al., 2019; Zhang et al., 2022). Therefore, we adopted PLS-SEM for further exploration.

When using PLS-SEM with normalized factors, multicollinearity diagnosis and model fit evaluation are essential for model quality evaluation. Multicollinearity diagnosis (Cenfetelli and Bassellier, 2009), which relies mainly on the variance inflation factor (VIF), is an important issue in model quality evaluation. High multicollinearity between factors may cause inaccurate path coefficient, leading to these factors being misinterpreted as unimportant or invalid (Hair et al., 2010). The model fit evaluation uses approximate fit indices like SRMR (standardized root mean square residual) and NFI (normed fit index). A desirable PLS-SEM is characterized by SRMR less than 0.08 and NFI greater than 0.90. To interpret the SEM results, a standardized path coefficient implies the direct effect of one factor on another, and its significance (p<0.05) can be evaluated by resampling procedures (Hair et al., 2011). The path coefficient was solved iteratively using ordinary least squares, and the significance (p<0.05) of each path coefficient was obtained by bootstrapping (5000 iterations). The indirect effect of one factor on another can be calculated through multiplying all indirect path coefficients, and the total effect is the sum of the indirect and direct effects.

4 Results

4.1 Accuracy assessment of the extracted snow phenology

Figure 3 shows the accuracy of the satellite-derived SP parameters compared with the ground-observed SP. The satellite-derived SCD, SOD, and SED values were in line with ground-observed ones, with R of 0.77, 0.95, and 0.97 (p<0.01), respectively. The biases of ground-observed and satellite-derived SCD, SOD, and SED were −1.75, −2.39, and 3.43 d, indicating that the satellite-derived SCD was less than the actual SCD, and the satellite-derived SOD was earlier than the actual SOD. In contrast, satellite-derived SED occurred later than the actual SED. Relevant studies have shown the bias of satellite-derived SP parameters ranges from 0 to 10 d (Chen et al., 2018; Hao et al., 2022; Wang et al., 2021). The biases (−2.39, 3.43, and −1.75 d) in our study were within this range.

Figure 3 The accuracy of SP parameters of (a) SCD, (b) SOD, and (c) SED evaluated by ground-observed values from 2002 to 2022. Note that SCD, SOD, and SED denote snow cover days, snow onset date, and snow end date, respectively. DOS represents the day of the snow season. ${}^{***}$ indicates significance at the level of 0.01.

4.2 Spatial patterns and temporal trends of snow phenology

The spatial patterns and temporal trends of SP at the pixel level are shown in Figs. 4 and 5, respectively. Figure 4a indicates the spatial distribution of the multiyear-averaged SCDs from 2002 to 2022, showing a highly uneven snow cover distribution across the plateau. Areas with an SCD of more than 60 d are considered stable snow cover areas and sources of water (Tang et al., 2022), accounting for approximately 17.07 % of the total TP and being primarily distributed in high-elevation mountain ranges. Some high-elevation areas had an SCD exceeding 180 d, accounting for approximately 7.40 % of the entire TP, distributed primarily in the western Kunlun, Himalaya, and Nyainqêntanglha mountains. Areas of low SCD (<20 d) accounted for 55.95 % of the total TP, mainly concentrated in low-elevation areas and the hinterland of the plateau. Overall, the pattern of SCD was in line with topography in a visual sense, displaying the regularity of a high SCD in high-elevation mountains and low SCD in the low-elevation plains. We also examined the temporal trend in SCD and its significance. As Fig. 5a shows, the SCD decreased in 50.59 % of the TP during the 20-year period and significantly decreased by 4.62 % of the total TP at a mean rate of −1.74 d yr⁻¹ (p<0.1), primarily in the west Kunlun, east-central Himalaya, and south Nyainqêntanglha mountains. In contrast, areas with significantly increased SCD accounted for 2.45 % of the total TP at a mean rate of 1.12 d yr⁻¹ (p<0.1) and were mainly scattered in the east-central part of the TP.

Figure 4 The spatial pattern of the multiyear-averaged SP of (a) SCD, (b) SOD, and (c) SED on the TP from 2002 to 2022. Note that pixels not detected for 20 years or identified as lake area were considered snow-free pixels (i.e., the white area in the figure).

Figure 5 Trends in SP of (a) SCD, (b) SOD, and (c) SED on the TP from 2002 to 2022. Note that the map at the bottom left of each panel indicates the significant level. The symbol “de” in the legend indicates a significantly decreased trend, and “in” indicates a significantly increased trend.

Distribution of SOD and SED over the TP also presented extreme spatial heterogeneity, which was visually consistent with elevation (Fig. 4b and c). DOS (day of the snow season) 1 is equivalent to 1 September of the previous year. SOD was mainly concentrated in DOS 90–180 (December–February), whereas SED was concentrated primarily in DOS 150–240 (February–April). In high-elevation mountain ranges, SOD appears earlier in general, whereas SED starts later. But in low-elevation areas such as the hinterland of the TP, SOD starts later and SED appears earlier. Not only are the variations in SOD and SED spatially complicated, but they also exhibited significant temporal heterogeneity (Fig. 5b and c). A total of 58.56 % of the total TP had delayed SOD, and 2.34 % of the TP had a statistically significant (p<0.1) delayed SOD at a mean rate of 2.90 d yr⁻¹ (Fig. 5b). In contrast, areas with significantly advanced SOD only accounted for 0.46 % of the total area at a mean rate of −3.69 d yr⁻¹ (p<0.1) and were sparsely distributed in the northwestern and eastern TP. For SED, 98.41 % of the TP showed no significant trend, 1.52 % showed a significant advanced trend at a rate of −2.49 d yr⁻¹ (p<0.1), and 0.07 % showed a significant delayed trend (p<0.1) (Fig. 5c).

4.3 Direct and indirect effects of important factors on snow phenology

Based on the related literature (Guo et al., 2022; Huang et al., 2020; Qi et al., 2021), 11 initial associated factors in four categories were selected: meteorological factors (air temperature, precipitation, wind speed, humidity, and shortwave radiation), geographical location (latitude and longitude), topography (elevation, slope, and aspect), and vegetation greenness (NDGI). First, multicollinearity diagnosis was performed on all initial factors, among which humidity did not pass the multicollinearity diagnosis (VIF = 10.8) and was excluded. Then, we constructed a SEM using the remaining factors and established multiple paths in the model. Non-significant and low coefficient paths were removed, and the total effect (TE) of each factor was calculated based on the final paths (Fig. 6). Factors with a small absolute TE should be removed because of their small contribution and limited explanations. An appropriate threshold was selected as the basis for factor elimination. Hair et al. (2011) indicated that there is no uniform standard for selecting thresholds in different research fields. Considering the rationality of factor selection and referring to the relevant literature (Shen et al., 2022; Zhang et al., 2022), we used 0.01 as the threshold for factor elimination. Aspect was eliminated for all SP (Fig. 6), and the final SEMs were constructed to discuss the direct and indirect effects of the remaining nine factors.

Figure 6 The absolute total effect of associated factors affecting SP of (a) SCD, (b) SOD, and (c) SED based on SEM.

Figure 7 presents the results of the final SEM and its standardized regression path coefficient (PC) (p<0.05). It indicates that there are complex interactions between SP and associated factors on the TP. Temperature, wind speed, NDGI, and shortwave radiation had negative effects on SCD and positive effects on SOD, while the other factors had positive effects on SCD and negative effects on SOD. Furthermore, temperature, wind speed, and shortwave radiation had negative effects on SED, while the remaining factors had positive effects on SED. Meteorological factors (air temperature, precipitation, wind speed, and shortwave radiation) exerted both direct and indirect effects on SP, whereas vegetation greenness (NDGI) mainly had a direct impact, and geographical location (latitude and longitude) and topography (elevation and slope) had an indirect effect on SP (Fig. 7). Meteorological factors were the dominant factors affecting SP, especially the temperature, which had the strongest effect on SCD, SOD, and SED, with the TE values of −0.447, 0.272, and −0.379, respectively. The influences of precipitation on SCD and SOD were also relatively significant (absolute value of TE >0.19). However, its impact on SED was limited (TE = 0.049). Wind speed exhibited a strong effect on SCD and SED (absolute value of TE >0.15), while its effect on SOD was relatively limited (TE = 0.062). The effect of NDGI on SCD and SOD was strongest (absolute value of TE >0.17). Geographical location and topographic conditions affected SP by determining meteorological conditions and the growth of vegetation. For example, elevation indirectly affected SCD by assessing the distribution of temperature (PC = −0.805), precipitation (PC = −0.261), shortwave radiation (PC = 0.330), wind speed (PC = 0.418), and NDGI (PC = −0.649); hence, the TE of elevation on the SCD was 0.290 (Fig. 7a). The influence of latitude on all SP parameters was greater than that of longitude.

Figure 7 The SEM of each SP parameter of (a) SCD, (b) SOD, and (c) SED. Note that the red line implies a positive effect, while the blue line denotes a negative effect. All path coefficients are statistically significant (p<0.05).

5 Discussion

Compared to traditional statistical methods, the SEM adopted in this study can accurately quantify the direct and indirect relationships between SP and associated factors. This makes it particularly well-suited for understanding the mediating effects of these factors. Based on the results of this model, the influences of relatively important factors are discussed below.

5.1 Response of snow phenology to meteorological factors

Meteorological factors (air temperature, precipitation, wind speed, and shortwave radiation) can directly affect SP as well as indirectly affect it by influencing the growth of vegetation, whereas the direct effect was much greater than the indirect effect. Among these factors, air temperature was the dominant meteorological factor. A temperature below 0 °C is conducive to the increase and maintenance of snow cover, resulting in longer SCD, earlier SOD, and later SED (Moran-Tejeda et al., 2013; Scalzitti et al., 2016). As another important factor, precipitation has been indicated by prior research to have a comparatively lower contribution to SP (Guo et al., 2022; Huang et al., 2020; Wang et al., 2021), which is inconsistent with the results of our research (Fig. 7). This discrepancy may be caused by the different resolutions of precipitation data. The spatial resolution of the precipitation data we used was approximately 3 km, while the resolution in previous studies was relatively rough (25 km or data from meteorological stations). The 3 km resolution precipitation data we used reveal greater details of precipitation distribution and provide more accurate precipitation information, thus obtaining a more reliable correlation between precipitation and SP.

To further explore the response mechanisms of SP to temperature and precipitation, their relationship was analyzed at the pixel level using Pearson's correlation coefficients (Fig. 8). SCD and SED exhibited a negative correlation with temperature across the majority of regions, accounting for 88.09 % and 73.81 % of the TP, respectively (Fig. 8a and e). Conversely, they demonstrated a positive correlation with precipitation in 59.67 % and 61.02 % of the areas of the TP (Fig. 8b and f). SOD was positively correlated with temperature in 78.02 % of the TP (Fig. 8c) and displayed a negative correlation with precipitation in 70.90 % of the TP (Fig. 8d). We also calculated the correlation coefficients of SP with temperature and precipitation for each elevation range (Fig. 8g, h, and i). The role of temperature and precipitation in SP varies with elevation. In areas with low elevation (<1500 m) and high elevation (>6700 m), the relationships between SP and temperature and precipitation were complex, with neither factor dominating significantly. For SCD, the influence of precipitation exceeded that of temperature and became more important in 1500–5400 m (Fig. 8g). However, temperature played a more crucial role in 5400–6700 m. Similarly, the effect of precipitation on SED was more significant in 1800–5400 m, after which the effect of temperature was more important in 5400–6700 m. The correlation between SOD and precipitation was stronger than its correlation with temperature across 4700–5500 m. Overall, we identified that the relative importance of temperature and precipitation shifts with elevation. Our findings differ from previous studies, where researchers generally found an elevation threshold where the importance of each changed. For instance, Moran-Tejeda et al. (2013) discovered a threshold elevation of approximately 1400 m, below which temperature was the primary explanatory variable of snow cover in Switzerland, and above which precipitation was a better predictor. Scalzitti et al. (2016) found a range of threshold elevations (1580–2181 m) that separated the importance of temperature from that of precipitation in the western United States. The reasons for such difference with our research could be as follows: (1) the complex terrain and different climate conditions of the TP and (2) the unique changes in temperature and precipitation at different elevations on the TP. The confounding influence of warming temperatures and varying precipitation remains a challenge in explaining the observed variations in the SP.

Figure 8 Correlation between temperature and SP of (a) SCD, (c) SOD, and (e) SED at pixel level; correlation between precipitation and SP of (b) SCD, (d) SOD, and (f) SED; and correlation coefficient between both temperature and precipitation and SP of (g) SCD, (h) SOD, and (i) SED under different elevation ranges.

This study still has some limitations in exploring the effect of meteorological factors on SP. We calculated the average (or total) metrics of each meteorological factor during the snow accumulation season to explore their relationship with SOD using SEM. Similarly, we calculated the average (or total) metrics of each meteorological factor during the snowmelt season to investigate their relationship with SED (Chen et al., 2018; Ma et al., 2023). The purpose of this is to address the time dependence of meteorological factors. In the future, we intend to explore more methods for effectively separating the time dependence of meteorological factors when studying their effects on snow phenology. Additionally, previous studies suggested potential lagging effects of meteorological factors on the SP of the TP. For instance, Li et al. (2019) found that anomalous precipitation could lead to subsequent snow cover variations on the TP with a delay of approximately 5 d. Ren and Liu (2018) observed a strong negative correlation between snow cover and the average temperature of 1 to 4 months prior to snow accumulation. Despite the potential importance of lagging effects on the snow cover of the TP, our SEM model requires concurrent observations of dependent and independence variables for accurate quantification of the causal relationships between them. Therefore, the lagging effects were neglected in this study.

5.2 Topography control on snow phenology

Previous studies have demonstrated the effect of complex topography on SP (Guo et al., 2022; Jain et al., 2008; Ma et al., 2020). Still, they did not reveal the specific action processes (direct or indirect) of various topographic factors affecting SP. This study indicated that topographic conditions indirectly affected SP by modulating meteorological conditions and the growth of vegetation using SEM. Elevation is the primary factor due to its vital effect on local microclimates, particularly in mountainous regions. In high-elevation mountains, the temperature and specific humidity are lower and snowfall is higher, creating favorable conditions for snowfall and maintenance. Owing to the strong blocking effect from large mountains, most areas of the hinterland of the TP have relatively scarce snow (Wang et al., 2017). Snow is more likely to slide on a steeper slope, which prevents snow accumulation, whereas a flatter slope is conducive to snow deposition. As shown in Fig. 7, slope had little effect on SOD (compared with SCD and SED), indicating that the snow accumulation was mainly influenced by climate and less by slope.

Elevation had a significant effect on SP, but whether the trends in SP from 2002 to 2022 on the TP was elevation-dependent was unknown. To investigate the elevation effect, we divided elevation into 50 m intervals, and the average trend values of statistically significant pixels (p<0.1) were calculated for each elevation category. Then, we selected elevation categories with more than 50 samples to analyze the correlation between trends in SP and elevation. Our finding reveals a strong negative correlation between trend in SCD and elevation, with a correlation coefficient of −0.73 (Fig. 9a). This negative correlation was strongest in 0–4000 and 4900–5900 m. A moderate positive correlation exists between the trend in SOD and elevation (R=0.59), and this positive correlation was most significant in 4100–5800 m (Fig. 9b). The correlation between SED and elevation was 0.33, which was not significant (Fig. 9c). Therefore, there exists a strong elevation dependence for the trend in SCD ($R=-\mathrm{0.73}$), a moderate elevation dependence for the trend in SOD (R=0.59), and no significant elevation dependence for the trend in SED (R=0.33) from 2002 to 2022. The elevation dependence of the SCD trend is consistent with the previous study of Ma et al. (2023). However, Ma et al. (2020) found no elevation dependence for both SOD and SED trends. This discrepancy is very likely because the meteorological data used in their study were mainly distributed in the eastern TP, while the western TP was much less considered. This further highlights the high spatial heterogeneity of SP in the TP.

Figure 9 Trends in SP of (a) SCD, (b) SOD, and (c) SED under different elevation ranges on the TP from 2002 to 2022 and the SEM based on delayed trend in SED and influencing factors (d).

The non-significant elevation dependence of SED is due to the competing effects of the delayed trend within 3800–4900 m elevation and the advanced trend of other elevation ranges (Fig. 9c). To explain the counterintuitive delayed trend in the context of regional warming, we further applied a structural equation model to explore the causal relationships. The results indicate that wind speed significantly influences the delayed trend in SED, exhibiting a strong positive correlation with a total effect of 0.497 (Fig. 9d). The delayed trend in SED becomes more pronounced as wind speed increases. This effect of wind speed on snow cover is mainly through the blowing snow process, which leads to increased snow accumulation through snow redistribution and thus a delayed trend in SED regionally. Li et al. (2012) also found a distinct occurrence of blowing snow at an elevation of 4146 m, where pronounced redistribution of snow cover occurred.

5.3 Other factors affecting snow phenology

The TP has a vast territory spanning multiple longitude and latitude zones, and longitude and latitude cannot be ignored in snow cover research. A greater longitude (i.e., further east) caused a shorter SCD, later SOD, and earlier SED (Fig. 4). From east to west, SCD decreased by 123 d, SOD was delayed by 77 d, and SED advanced by 48 d. The elevation of the TP decreases from west to east, resulting in a correlation between elevation and longitude (R=0.532; Fig. 10a). This correlation can cause the effect of longitude on SP folded under the effect of elevation on SP. However, controlling the influence of elevation, we still find the longitudinal dependency of snow cover. Taking the SCD as an example, at a fixed elevation (e.g., 5000 m), a correlation still exists between longitude and SCD (R=0.356; Fig. 10b). This implies that longitude could be interacting with additional factors to shape the spatial and temporal distribution of snow cover on the plateau. Vegetation greenness can change snow redistribution by intercepting the snow and regulating the balance of the solar radiation budget (Domine et al., 2016), and vice versa, the processes and patterns of snow accumulation and melting affect vegetation greenness, species distribution, and community structure (Barrere et al., 2018). Most previous studies explored the response of vegetation greenness to snow cover, and the role of vegetation greenness on snow cover has rarely been investigated (Qi et al., 2021; Xu et al., 2022b). Our study showed that vegetation had mainly direct effects on snow cover, with particularly pronounced effects on SCD and SOD.

Figure 10 Scatter plot of (a) elevation and longitude and (b) longitude and snow cover days (at a fixed elevation of 5000 m).

The dynamics of snow accumulation and melting are influenced by various factors. In addition to the factors analyzed in this paper, other factors may also play important roles in SP. Ground temperature primarily influences the structure and stability of snowpack by regulating energy exchange at the soil–snow interface (Rixen et al., 2022). As the ground temperature increases, the substrate absorbs additional thermal energy, which is conveyed to the base of the accumulated snow through heat conduction, resulting in melting of the lower snow layers. The heat from the warmer soil in the snow-free area can be transferred to the colder soil below the snow-covered area. Liquid water can also be transferred from the snow-free soil to the snow-covered soil, thus melting snow (Fassnacht et al., 2006). Therefore, soil properties (e.g., soil moisture) can also affect snow cover. In addition, atmospheric pollutants, especially those referred to as light-absorbing aerosols, such as black carbon, brown carbon, and dust, can warm the atmosphere (Kang et al., 2019; Ji et al., 2015). After being deposited onto snowpack, these light-absorbing particles can reduce the surface albedos of snowpack and promote its melting (Zhang et al., 2018; Lau and Kim, 2018). Despite the potential importance of these factors to the SP of the TP, they were not analyzed in this study due to limited data availability of these factors over extended spatial and temporal scales. Developing continuous high-resolution, spatiotemporal datasets for these factors will be useful in future efforts to comprehensively quantify the response of SP to changing climate conditions.

6 Conclusions

In this study, we investigated the spatiotemporal variability in SP on the TP from 2002 to 2022 based on the daily gap-free HMRFS-TP dataset. We also quantified the direct and indirect effects of the associated factors on SP. The spatial patterns of SP were vertically zonal, with higher elevations having a longer SCD, earlier SOD, and later SED. SCD had a statistically significant (p<0.1) decreased trend in 4.62 % of the TP at a mean rate of −1.74 d yr⁻¹. SOD was significantly delayed in 2.34 % of the TP at a mean rate of 2.90 d yr⁻¹ (p<0.1), while SED was significantly advanced in 1.52 % of the TP at a rate of −2.49 d yr⁻¹ (p<0.1). Additionally, there exists a strong elevation dependence for the trend in SCD ($R=-\mathrm{0.73}$), a moderate elevation dependence for the trend in SOD (R=0.59), and no significant elevation dependence for the trend in SED (R=0.33). Meteorological factors can directly affect SP as well as indirectly affect it by influencing the growth of vegetation, whereas the direct effect was much greater than the indirect effect. Geographical location and topographic conditions indirectly affected SP by modulating meteorological conditions and the growth of vegetation. Vegetation primarily influences SP by intercepting the snow and regulating the balance of the solar radiation budget. As two rather important factors, we identified that neither temperature nor precipitation showed consistently high importance as elevation increased.

This study explored the dynamic variation in snow cover and revealed the mediating effects of multiple factors in its changing process, which contributed to providing a strategic basis for predicting and solving the problems of climate change, hydrological cycle, and ecological balance in the future in the context of global warming on the TP.