*3.2. Methods*

### 3.2.1. SD Data Downscaling

We define the spatial resolution of 0.25◦ as the lower resolution and that of 500 m as the finer resolution. The purpose of this is to downscale the lower resolution SD data to the finer resolution by considering and integrating multiple influential factors.

(1) Statistical downscaling

Pilot correlation analysis revealed that topographic and geographic factors such as elevation, slope, aspect, longitude, and latitude have significant effects on SD; they affect or regulate climatic process such as snowfall at local or even larger scales [54,55]. All the above factors were calculated and resampled into the lower and finer spatial resolutions. Taking the influential factors at coarse resolution as environmental variables, a multiple linear regression method is used to calibrate the statistical formula based on the original monthly 0.25◦ SD data (*SDa*) as the target variable. Application of the formula results in a simulated series of SD at both the lower (*SDb*, Equation (1)) and finer (*SDc*, Equation (3)) spatial resolutions. Residuals (Equation (2)) of *SDb* to SD*a* (Δ*SDlow*) are resampled into the finer 500 m resolution (Δ*SDhigh*) using bilinear interpolation. The sums of *SDc* and Δ*SDhigh* are the statistically downscaled SD data (*SDd*, Equation (4)).

$$SD\_b = A \times X\_{1,low} + B \times X\_{2,low} + C \times X\_{3,low} + D \times X\_{4,low} + E \times X\_{5,low} + F \tag{1}$$

$$
\Delta SD\_{low} = SD\_a - SD\_b \tag{2}
$$

$$SD\_c = A \times X\_{1,high} + B \times X\_{2,high} + C \times X\_{3,high} + D \times X\_{4,high} + E \times X\_{5,high} + F \tag{3}$$

$$SD\_d = SD\_\emptyset + \Delta SD\_{\text{high}} \tag{4}$$

where *Xn* (*n* = 1, 2, 3, 4, 5) represent raster data of elevation, slope, aspect, longitude, and latitude, respectively. Subscripts *low* and *high* denote the lower and finer resolutions in space, respectively.

(2) Improvement based on snow cover data

The original resolution of the SD data is 0.25◦, and although it has been geographically and topographically corrected through statistical downscaling, the inaccurate influence is still present. For example, SD occurs in low-altitude regions in warmer months, a situation that seldom occurs, as common sense in the study area dictates. By introducing the spatial distribution probability function of snow, the 500 m MODIS snow cover data are used to modify and improve the precision of the statistical downscaling. Given values of 1 and 0 representing snowy and snowless pixels, respectively, the original MODIS snow cover raster data are converted into binary ones. We define the period from September 1 of one year to August 31 of the next year as a snow hydrological year (SHY), and the accumulation days of snow (*ADS*) in each SHY are calculated from September 2002 to August 2018 at a spatial resolution of 500 m. Then, the cumulative days of snow (*CDS*) in a domain containing 55 × 55 *ADS* grids are counted, approximately corresponding to the spatial resolution of the original SD data (*SDa*). The snow distribution probability (*P*) can be determined as follows:

$$P = \frac{ADS\*3025}{CDS} \tag{5}$$

The product of *SDd* and *P* (Equation (6)) is what objective (1) aims at, which is an improvement of the statistically downscaled SD based on snow cover data, comprehensively reflecting a high-resolution SD (*SDfin*) controlled or influenced by geography, topography, and snow distribution (Figure 2).

$$SD\_{fin} = SD\_d \* P \tag{6}$$

**Figure 2.** Schematic diagram of the statistical-based method and the improvement from introducing snow distribution probability for downscaling SD data.

### 3.2.2. LSSP Indicator Extraction

Snow cover maintaining days (SCD), snow cover start date (SCS), and snow cover melt end date (SCM) are calculated based on the binary retreated MODIS snow cover data. Among them, SCS and SCM are important LSSPs determining SCD [56,57], which represent the dates when a monitored pixel starts accumulating snow and ends melting in a SHY. SCD is the number of days that each pixel is covered by snow in a SHY. The larger the SCD is, the longer and the more the snow reserves.

$$\text{SCS} = F\_d - \text{SCD}\_{bFd} \tag{7}$$

$$\text{SCM} = \text{SCD}\_{aFd} + F\_d \tag{8}$$

$$\text{SCD} = \sum\_{i=0}^{N} (\mathcal{S}\_i) \tag{9}$$

In the above equations, *Fd* is a fixed number representing the date when the largest snow cover occurs during the period from 2002 to 2018. The statistics resulted in the *Fd* date being 12 January 2008, the 134th day in the SHY. *SCDbFd* and *SCDaFd* represent the numbers of snow cover days before and after *Fd*, respectively, in each SHY. *N* is the upper limit for a specified time range, valued as 1 to represent a complete SHY; *Si* is the binary retreated pixel value of daily snow cover (snowy or snowless).

### 3.2.3. Trend Analysis

Sen's slope is selected for series variation amplitude statistics [58,59]. The method can reduce or prevent the impact of data anomalies and omissions when evaluating the trend and range of time series changes. An orderly column is constructed with the change rates of sample sequences of different lengths. Variable testing is then performed statistically according to the given significance level to obtain the value range of the change rates, and the median is used to determine the variation trend and magnitude. The equation is as follows: 

$$SS\_{i\bar{j}} = \text{MEDIAN} \frac{(X\_{\bar{j}} - X\_{i})}{(j - i)} \tag{10}$$

where *SSij* is Sen's slope, *Xi* and *Xj* represent the sequential values corresponding to time *i* and *j*, respectively, where 1 < *i* < *j* < *n* and *n* is the length of the series.

The Mann–Kendall method [60,61] is a nonparametric test approach used to determine the significance of a trend analysis [62].

$$S = \sum\_{i=1}^{n-1} \sum\_{j=i+1}^{n} \begin{cases} 1 & y\_i - y\_j > 0 \\ 0 & y\_i - y\_j = 0 \\ -1 & y\_i - y\_j < 0 \end{cases} \tag{11}$$

$$Z = \begin{cases} \frac{S - 1}{\sqrt{s(S)}} & S > 0 \\ 0 & S = 0 \\ \frac{S + 1}{\sqrt{s(S)}} & S < 0 \end{cases} \tag{12}$$

where *yi* and *yj* represent the snow phenology indicators in SHY *i* and *j*, respectively; *n* represents the length of the sequence. A positive value of the statistic *S* indicates an increasing trend of the data series, while a negative one indicates a decreasing trend of the series. The value of *Z* is in the range of ( − <sup>∞</sup>, + ∞); for a given confidence interval *α*, if | *Z*| ≥ *Z*1−*<sup>α</sup>*/2, it indicates that there is a significant trend in the data series at confidence level α.

### 3.2.4. Analysis of Climate-Driven Influences

Correlation-based significance statistics, represented by partial correlation coefficients, are used to examine the strength of climate influences on regional LSSP. To be specific, assuming there are *j* (*j* > 2) variables (*Xl*1 , *Xl*2 ...... *Xlg* , *Xlh* , *Xli* , *Xlj*), the formula for the (*j* − 2)th order partial correlation coefficient of any two variables *Xli*and *Xlj*is:

$$r\_{l\_i l\_j \cdot l\_1 l\_2 \cdots l\_h} = \frac{r\_{l\_i l\_j \cdot l\_1 l\_2 \cdots l\_{h-1}} - r\_{l\_i l\_h \cdot l\_1 l\_2 \cdots l\_{h-1}} r\_{l\_j l\_h \cdot l\_1 l\_2 \cdots l\_{h-1}}}{\sqrt{\left(1 - r\_{l\_i l\_h \cdot l\_1 l\_2 \cdots l\_{h-1}}^2\right) \left(1 - r\_{l\_j l\_h \cdot l\_1 l\_2 \cdots l\_{h-1}}^2\right)}} \tag{13}$$

where *r* denotes the correlation coefficient, factors on the right side of the formula are the (*j* − 3)th order partial factors.

Four key climatic factors, including precipitation (P), land surface net solar radiation (SSR), maximum air temperature (T*max*), and minimum air temperature (T*min*), are selected to investigate the climate-driven strength of the LSSP, the 3rd order partial correlation coefficient is adopted. A *t* test is used to analyze whether the partial correlation coefficient between LSSP and climatic factors pass the 0.01 significance level (Table 1). The calculation formula is:

$$\mathbf{t} = \frac{r\_{l\_4l\_5 \cdot l\_1l\_2l\_3}^2}{1 - r\_{l\_4l\_5 \cdot l\_1l\_2l\_3}^2} \times \frac{n - m - 1}{m} \tag{14}$$

where *rl*4*l*5·*l*1*l*2*l*<sup>3</sup> is the partial correlation coefficient, *m* is the number of independent variables, and *n* is the sequence length. The partitioning criteria are listed in Table 1.



Note: *t* test-P, *t* test-SSR, *t* test-T*max,* and *t* test-T*min* represent the t significance test of LSSP with P, SSR, T*max*, and T*min*, respectively. [P] and [SSR] indicate that LSSP is driven by P or SSR. [P + SSR] indicates that LSSP is driven by both P and SSR. [P + SSR + T*max* + T*min*] means that LSSP is conjointly driven by P, SSR, T*max*, and T*min*. [NC] means nonclimate-driven. There are four key climatic factors and a total of 16 driving force combinations or single factors. Not all of them are listed due to space limitations.

### 3.2.5. Sensitivity Analysis

The response of snow variation to climate change is diagnosed using the sensitivity coefficient [63]. The method is widely used in contribution separations of influential factors on hydrological processes. The sensitivity coefficient is calculated as:

$$\varepsilon\_x = \frac{\overline{\overline{x}}}{\overline{\overline{y}}} \times \frac{\sum\_{i=1}^{n} (x - \overline{x})(y - \overline{y})}{\sum\_{i=1}^{n} (x - \overline{x})^2} \tag{15}$$

where *εx* is the sensitivity coefficient of *y* (LSSP) to *x* (climatic factors), indicating that the *<sup>ε</sup>x*% change in LSSP is caused by 1% variation in a climatic element. *x* and *y* are the multiyear averaged values of *x* and *y*. In the following statement, *εa*−*b* represents the sensitivity of LSSP to climatic factors, *a* is climatic factors such as P, SSR, T*max*, and T*min*, *b* is LSSP indicators such as SD and SCD.

## **4. Results**

### *4.1. Evaluation of the Improved SD Downscaling Method*

Gauge-based SD observations at 10 stations in and near the SGP (Figure 1c) were adopted as references, based on which absolute errors were calculated using SD data before and after downscaling, as shown in Figure 3. The vertical ranges cover all errors, including both the positive and negative ones, while transverse widths indicate occurrent frequencies. It can be seen that the positive downscaled SD (Δ*SDfin*) errors at all stations

and the negative errors at most of the 10 stations are reduced, indicating an effective optimization for the elimination of both the over- and underestimations of the initial SD data. Frequencies of SD error valued at 0 were found to be the highest at all stations, and differences were not clear between the two sets of data, although an overall but slight decrease appeared. As a whole, using the improved downscaling method, the spatial resolution and real representation of SD were verified to be acceptable, and the downscaled data were satisfactory for LSSP analysis in the SGP region.

**Figure 3.** Evaluation of the improved downscaling method based on the differences from the gaugebased SD observations. Δ*SDa* and Δ*SDfin* represent the absolute error between the SD data before and after downscaling, respectively. Subplots (**<sup>a</sup>**−**j**) correspond to the gauge stations for snow observation. Names of the stations labeled in the right corners of the plots.

### *4.2. LSSP Characteristics*

The binary retreated MODIS snow cover data were used to calculate time-based LSSP values, including the SCS, SCM, and SCD. Together with the downscaled SD data, a spatiotemporal analysis of LSSPs of the SGP from 2002 to 2018 was conducted.

### 4.2.1. Spatiotemporal Distributions

LSSPs present remarkable heterogeneity in time and space on the SGP. High SD values are mainly found in areas near or at mountain divides, while low SD values are generally distributed in valley areas featuring lower altitudes. For example, in the southwestern and southeastern parts where the main stream of the Yellow River and the Bailong River develop, SD can exceed 18 cm in high altitude zones but less than 10 cm in valleys (Figure 4a). The SCS and SCM show opposite distribution patterns in value space (Figure 4b,c). In or near mountain divides, snow begins to accumulate early (i.e., the earliest is the 34th day in a SHY) but ends melting late (i.e., after the 255th day in a SHY). In most areas, the SCS starts after the 120th day, corresponding to SCM days before the 195th day, and both result in a relatively short SCD time in a SHY, i.e., most of the SCD is quantified to be less than 15 days in the SGP. Overall, the SD, SCS, SCM, and SCD spatially match well in phenology. For example, in areas such as mountain divides featuring deeper snow, accumulation starts earlier, melting ends later, and snow cover lasts longer.

The SD is larger in months from November to the following March than that in the rest of a SHY, the maximum value generally occurs in February, with a multiyear average value of 4.03 cm (Figure 5a). For the presence of snow cover, SCD presents a "multipeak" pattern in a SHY (Figure 5b). Generally, snow cover lasts the longest in November of one year and March of the following year, especially in March, when the existence of snow can reach 8 days. In the context of the regional climate, there is rarely snow cover in the SGP from June to August.

**Figure 4.** The spatial distribution of LSSPs on the SGP. Subplots (**<sup>a</sup>**–**d**) are for SD, SCS, SCM, and SCD.

**Figure 5.** Multiyear averaged monthly distribution of SD (**a**) and SCD (**b**) in SGP (2003–2018).

4.2.2. Spatiotemporal Variations

LSSP variation was found to interannually fluctuate according to quantification of the representative index, including SD, SCS, SCM, and SCD. Generally, higher values of SD and SCD occurred simultaneously, corresponding to smaller SCS values and larger SCM values (Figure 6). During the period from 2003 to 2018, the value ranges of the LSSP index, such as the SD, SCS, SCM, and SCD, were quantified from 2.67 to 10.53 cm, 113 to 134 d, 142 to 188 d, and 16 to 75 d, averaged to 5.15 cm, 123 d, 158 d, and 34 d, respectively. During

the statistical period, the maximum SD was found in 2009, when the SCS and SCM were 121 d and 181 d in the SHY, respectively, corresponding to an SCD length of 60 d.

**Figure 6.** Interannual variations of the four selected LSSP indices on the SGP (2003–2018).

The variation trends and amplitudes of the LSSP during the time period from 2003 to 2018 were measured spatially on the SGP based on Sen's slope method, as shown in Figure 7. A decreasing trend of SD in most of the area occurred, especially in the southwestern part, which is located in the main channel of the Yellow River. A relatively small area had an increase in SD, as in the southeastern part, which dominates the upstream section of the Bailong River. The variation in SD was not significant at *p* ≥ 0.01 on most of the SGP, especially in the central and northern parts, similar to the significance test for the spatial statistics of the other three LSSP indices. The areal percentage of the reduced SD occupied 60.36%, and the variation amplitude was overall quantified into a rate of −0.06 cm/a across the whole SGP (Figure 7a). Areas with high altitudes showed significant changes in the SCS and SCM, such as in or near mountain divides. In terms of the areal percentage, areas with increased SCS accounted for 75.10%, indicating that the start date of snow accumulation on most of the SGP presented a delay. Areas with decreased SCM accounted for 65.78%, indicating that most of the snow melt ended earlier. However, the aggregative variation was rated as 0.53 d/a and 0.45 d/a for SCS and SCM, respectively (Figure 7b,c). The amplitude of the former was greater than that of the latter, resulting in an overall reduction in SCD of −0.37 d/a, which was specifically significant in the southwestern and southeastern parts of the SGP. The statistics resulted in a larger areal percentage of 56.90% for the increase in SCD, and the overall reduction was due to the magnitude of the decrease (Figure 7d).

### *4.3. Terrain Influence on LSSP*

Terrain factors, including altitude, slope, and aspect, were calculated and reclassified using value intervals of 500 m, 10◦, and 45◦ for analysis of the topographic influence on the LSSP. Altitude zones between 3000~3500 m and 3500~4000 m are mostly distributed in the region (Figure 8a), while few of the slopes are more than 50◦ (Figure 8b). The aspect differentiation is not significant, with relatively more differentiation in the northeast and less differentiation in the southeast and south (Figure 8c).

**Figure 7.** Spatial variation of the four selected LSSP indices across the SGP (2003–2018). Subplots (**<sup>a</sup>**–**d**) are for SD, SCS, SCM, and SCD in order). "∗" indicates an area with significantly tested LSSP series (*p* ≤ 0.01). Fan diagrams in the upper left corners of the subplots indicate areal proportions of increasing and decreasing, differentiated using blue and orange, while regional averages are indicated by the white part in the middle.

### 4.3.1. Terrain-Based Distribution

LSSP presented a similar value distribution corresponding to altitude and slope. Statistically, the higher the altitude, or the steeper the slope, the greater the SD (Figure 9a,c), SCM (Figure 9c,f) and SCD (Figure 9d,h), and the smaller the SCS (Figure 9b,g). In particular, the variation amplitudes of the four LSSP indices increased along with the above two terrain factors, indicating that LSSP variability strengthened in regions with higher altitudes and steeper slopes. Generally, snow does not easily accumulate on steep slopes, our analysis showed an inverse pattern. It was found that steep slopes are mostly located in regions with high altitudes on the SGP and are more prone to snowfall, which might be the reason why the higher SD values were more distributed [64]. The influence of aspects on LSSP were statistically close due to the relative equilibrium at the regional scale. South-facing areas, including the southeast, south, and southwest facings, had smaller SD, SCM, and SCD, and larger SCS (Figure 9i–l), reflecting a lower probability of snow accumulation on sunny slopes.

**Figure 8.** Reclassification of the three selected topographic factors together with their areal statistics across the SGP. Subplots (**<sup>a</sup>**–**<sup>c</sup>**) are for altitude, slope, and aspect, respectively.

**Figure 9.** Statistics of LSSP based on terrain. Subplots (**<sup>a</sup>**–**d**), (**<sup>e</sup>**–**h**) and (**i**–**l**) correspond to altitude, slope, and aspect, respectively.

### 4.3.2. Terrain-Based Variation

Statistics were conducted to diagnose the variation difference in the LSSP corresponding to the reclassified terrain. The results revealed that the decrease in SD mainly occurred in regions with altitudes higher than 3000 m, while in other altitude zones, the SD increased. In zones with slopes ranging between 70~80◦ or less than 40◦, the SD was found to decrease, although an increase occurred in other slope zones. All aspects, correspondingly exhibiting almost all altitudes and slopes in zones, were found to decrease with SD. In terms of the time phenology, an increase in SCS, together with a decrease in SCM, was found in most of the altitude zones, resulting in a general decrease in SCD or a shortened snow-maintaining time, especially in regions with altitudes higher than 3000 m, consistent with the variation in SD (Figure 10a). The SCS in all slope zones presented a delay (increase), while the SCM presented a decrease (earlier), except in zones with slopes less than 20◦ (Figure 10b). Statistics based on aspects resulted in an overall increase in the SCS and SCM, although the variation rates of the latter were less than those of the former (Figure 10c). According to the above, given an altitude of 3000 m and a slope of 20◦, LSSPs on the SGP generally presented a pattern of "decrease in higher and steeper areas corresponding to a shortened snow duration, increase in lower and flatter areas corresponding to a relatively lengthened snow duration", which might have a profound relationship with the regional climate change under the warming background.


**Figure 10.** Quantified annual variation in SD (cm/a), SCS (d/a), SCM (d/a), and SCD (d/a) corresponding to the terrain factors of the SGP. Subplots (**<sup>a</sup>**–**<sup>c</sup>**) represent zonal statistics based on the reclassified altitude, slope and aspect, respectively.

### *4.4. Climate Influence on LSSP*

Analysis revealed that the SGP experienced an overall reduction in snow accumulation during the time period from 2003 to 2018, the duration of snow cover moved backward in a SHY. The above, represented by the SD and SCD indices (determined by SCS and SCM), can be considered comprehensive reflections of snow dynamics. Reasons leading to snow variation vary, in which the climate influence plays an important role [25,65]. In particular, regional precipitation (P), surface net solar radiation (SSR), maximum air temperature (T*max*), and minimum air temperature (T*min*) are key factors [66] and are thus selected to analyze the responses of the LSSP to regional climate change on the SGP.

### 4.4.1. Partial Correlation Analysis

The influence of climate change on regional LSSP is mechanically complex. To be more focused, the multicollinearity analysis between factors of climate and terrain are ignored at this stage because it is out of the present study's scope. Given the four selected climatic factors, the singular or synergistic effectiveness of multiple factor combinations on LSSP variation needs to be clarified. Coupling partial correlation (Equation (14)) and significance analysis (Table 1) can help. Among the 16 factor/combinations, the top five categories satisfying the t test standard were sorted by significantly influenced area. The SD order list was [SSR, T*min*, T*max*, SSR + T*min*, P] (Figure 11a), while that of the SCD was [T*max*, SSR, P, P + SSR + T*max*,P+T*max*] (Figure 11b), indicating different impacts of climate on regional LSSP variation, e.g., SD varied remarkably by SSR and T*min* effectiveness, although SCD was more susceptible to synergistic influences, mechanically showing more

complexity. Effectiveness of climatic driving factors on SD and SCD were found to be different, which may be due to local terrain regulations to different LSSP. For example, in the northern part of the SGP featuring lower altitude, the SD is more significantly related to lower air temperature that facilitating more snowfall, while the occurrence of factors with influence that can significantly drive the SCD is less because snow is usually sparse and shortly maintained there. In the southwestern part, the altitude is relatively higher and the air temperature is always lower; the abundant SSR, due to the high altitude, influences SD more, although SCD is dramatically affected by the higher air temperature because the snow event occurs commonly, especially in cold seasons. In summary, the singular effectiveness of the factors was verified to influence the LSSP more over a larger area, which may be related to the zonal differentiation of terrain and climate on the SGP. In other words, key factors affecting regional LSSP might tend to be single, which is in good agreemen<sup>t</sup> with the following sensitivity analysis.

**Figure 11.** Significance-based diagnosis of the relationship between LSSP and key climatic factors and their combinations on the SGP. Subplots (**<sup>a</sup>**,**b**) correspond to SD and SCD, respectively. Histograms in the left corners show the top five areas influenced (km2) and NC indicates an area that failed the significance test.

### 4.4.2. Sensitivity Analysis

Sensitivity (coefficient) analysis quantifies the percentage variation in the target variable (SD/SCD) caused by a 1% change in the environmental variable (the climatic factors). *εP*−*SD*, *εSSR*−*SD*, *εTmax*−*SD* and *εTmin*−*SD* represent the sensitivity of SD to P, SSR, T*max*, and T*min*, respectively. The results revealed that SD positively or negatively responded to all four climatic factors. Positive facilitation of SD due to increased P was widely characterized on the SGP, accounting for 67.74% of the area (Figure 12a), while the other three factors produced damping effects on SD; along with the increase in SD, the proportions of the area with negative impacts were above 80% (Figure 12b–d). The regional averages of *εP*−*SD*, *εSSR*−*SD*, *εTmax*−*SD* and *εTmin*−*SD* were 0.71%, −4.12%, −3.22%, and −2.95%, respectively, corresponding to SD variation rates when the climatic factor increased by 1%. Comparatively, SD was more sensitive to SSR on the SGP.

**Figure 12.** Spatial distribution of the SD sensitivity to the four climatic factors. Subplots (**<sup>a</sup>**–**d**) correspond to P, SSR, T*max* and T*min*, respectively. Fan diagrams in the upper left corners of subplots indicate the areal proportion of positively or negatively influenced and are differentiated using blue and orange, while the regional averages are the white part in the middle.

*εP*−*SCD*, *εSSR*−*SCD*, *εTmax*−*SCD*, and *εTmin*−*SCD* represent the sensitivity of SCD to P, SSR, Tmax, and Tmin, respectively. Similar results were obtained by calculations (Figure 13). In contrast, the regional averages of *εP*−*SCD*, *εSSR*−*SCD*, *εTmax*−*SCD*, and *εTmin*−*SCD* were 1.16%, −3.69%, −3.45%, and −4.01%, respectively, corresponding to SCD variation rates when the climatic factor increased by 1%. Comparatively, SCD was more sensitive to air temperature, especially the lower ones.

Overall, the LSSP on the SGP was influenced by and showed more sensitivity to solar radiation and air temperature. Local precipitation, especially snowfall, provides the matter base for snow formation, distribution, and accumulation. The low sensitivity might be related to the spatiotemporal correspondence and the interactive adaptation between terrain and climate.

### *4.5. Integrated Effectiveness of Climate and Terrain*

To explore the integrated effectiveness of climate and terrain, statistics of LSSP sensitivity (coefficients) to the climate were conducted based on the terrain. The results showed negative values of *εP*−*SD* with respect to the area, mainly at altitudes lower than 3000 m, while positive values were observed in the others, indicating that precipitation tended to positively contribute to SD in higher altitude regions [67]. Values of *εP*−*SCD* were found to be positive in most of the study area except in the altitude range of 1500~2000 m, indicating that it is not easy for snow to form and be maintained in lower altitude regions. The sensitivity coefficients of the two LSSP indices to SSR, T*max*, and T*min* were all negative, indicating that an increase in the three climate factors led to a decrease in snow in both

amount and persistence at all altitudes (Figure 14a,b). In areas with slopes greater than 30◦, the SD sensitivity to precipitation was found to be negative, indicating that snow hardly accumulated on steep slopes. The SCD sensitivity to precipitation was positive but decreased with increasing slope, indicating that the steeper the slope is, the smaller the effect of precipitation on snow persistence. The absolute value of LSSP sensitivity to SSR and T*max* decreased with increasing slope, while the trend of LSSP sensitivity to T*min* was opposite, meaning that snow on steeper land is less sensitive to solar radiation and higher air temperature, but more sensitive to lower air temperature (Figure 14c,d). LSSP sensitivities to P and SSR were not found to be significantly different in various aspects, while LSSP sensitivity to T*max* was found to be stronger on sunny (southward) slopes and LSSP sensitivity to T*min* was stronger on shaded (northward) slopes (Figure 14e,f).

**Figure 13.** Spatial distribution of SCD sensitivity to the four climatic factors. Subplots (**<sup>a</sup>**–**d**) correspond to P, SSR, T*max* and T*min*, respectively. Fan diagrams in the upper left corners of subplots indicate the areal proportion of positively or negatively influenced and are differentiated using blue and orange, while the regional averages are indicated by the white part in the middle.

**Figure 14.** LSSP sensitivity statistics based on combinations of climatic and terrain factors. Subplots (**<sup>a</sup>**,**b**) are altitude-based, (**<sup>c</sup>**,**d**) are slope-based, and (**<sup>e</sup>**,**f**) are aspect-based.

### **5. Discussion**

Located in a midlatitude inland region featuring limitations of circulation transportation, precipitation, especially snowfall in the cold season on the SGP, is relatively scarce [41]. Combined with the influence of land surface topography, snow cover presents remarkable fragmentation [68,69], leading to significant heterogeneity in LSSP variation [70]. For a better understanding, data with a fine spatial resolution are thus essential. To address this issue, we regressively analyzed the relationship between SD and the geographic and topographic factors on the SGP to statistically downscale the SD data product from a coarse spatial resolution of 0.25◦ to a finer one of 500 m. To eliminate the distortion of the original data, we introduced the MODIS-based spatial probability of snow cover for correction [7,71]. The coupling of the two ideas not only helped advance the downscaling method but also improved the data accuracy of snow depth, which can be considered an example to synthesize methods in data treatment.

Snow formation, accumulation, and redistribution are closely related to land surface terrain [72,73]. We found that on the SGP, the higher the altitude is, the larger the LSSP variability, indicating that snow is more significantly affected by environmental factors at higher altitudes [74]. The slope, together with aspect, is profoundly connected with incidence angle and eolian activities, acting on the amount of solar radiation and energy flux, influencing land surface energy budget, and affecting snow accumulation and melt rates [68,75]. Our study demonstrated the difference in LSSP in terrains with high complexity.

It was found that precipitation and air temperature are the key factors that most influence the regional LSSP [76]. As the basic material source for forming snow, precipitation cannot have an effect on the LSSP unless the air temperature is below the freezing point [77,78], which might be the explanation for the lower sensitivity of the LSSP to precipitation and the more remarkable sensitivity to air temperature [79,80]. Generally, absorption of solar radiation by snow cover relates to land surface albedo [81,82], during which the phase variation of the solar altitude angle in winter and spring plays an important role, dramatically affecting snow melt in a year [83,84]. Our study of the SGP agree with the above findings. The influence of key climatic factors on LSSP was further distinguished, and considerable differences in LSSP in response to the maximum and minimum air temperatures in time and space were quantitatively presented.

In terms of LSSP variation, the chemical composition and physical state of the atmosphere (e.g., aerosol load), structure and density of snow (e.g., crystal size) [85], variation in

environmental conditions (e.g., permafrost degradation [82]; land use/cover change [86]), artificial disturbance [87,88], and so on, directly or indirectly have profound impact on snow accumulation and melting and may cause LSSP variability [89]. In particular, a long period with snow is generally influenced by combinations of interactively functioning factors [90]. Our study examined individual effectiveness based on sensitivity analysis, and coupling studies were not conducted. Additionally, LSSP does not vary synchronously with climate, and delays with time often occur [78,91]. Furthermore, land surface vegetation also influences the LSSP [92]. The above processes, together with their effects on regional LSSP, have been neglected and may lead to uncertainties in the present study. It is hoped that shortcomings will be overcome in the future based on the improvement of data and methods.
