Trends and variability in snowmelt in China under climate change

Snowmelt is a major fresh water resource, and quantifying snowmelt and its variability under climate change is necessary for the planning and management of water resources. Spatiotemporal changes in snow properties in China have drawn wide attention in recent decades; however, country-wide assessments of snowmelt are lacking. Using precipitation and temperature data with a high spatial resolution (0.5^′; approximately 1 km), this study calculated the monthly snowmelt in China for the 1951–2017 period, using a simple temperature index model, and the model outputs were validated using snowfall, snow depth, snow cover extent and snow water equivalent. Precipitation and temperature scenarios developed from five CMIP5 models were used to predict future snowmelt in China under three different representative concentration pathway (RCP) scenarios (RCP2.6, RCP4.5 and RCP8.5). The results show that the mean annual snowmelt in China from 1951 to 2017 is 2.41×10¹¹ m³ yr⁻¹. The mean annual snowmelt values in Northern Xinjiang, Northeast China and the Tibetan Plateau – China's three main stable snow cover regions – are 0.18×10¹¹, 0.42×10¹¹ and 1.15×10¹¹ m³ yr⁻¹, respectively. From 1951 to 2017, the snowmelt increased significantly in the Tibetan Plateau and decreased significantly in northern, central and southeastern China. In the whole of China, there was a decreasing trend in snowmelt, but this was not statistically significant. The mean annual snowmelt runoff ratios are generally more than 10 % in almost all third-level basins in West China, more than 5 % in third-level basins in North and Northeast China and less than 2 % in third-level basins in South China. From 1951 to 2017, the annual snowmelt runoff ratios decreased in most third-level basins in China. Under RCP2.6, RCP4.5 and RCP8.5, the projected snowmelt in China in the near future (2011–2040; mid-future –2041–2070; far future – 2071–2099) may decrease by 10.4 % (15.8 %; 13.9 %), 12.0 % (17.9 %; 21.1 %) and 11.7 % (24.8 %; 36.5 %) compared to the reference period (1981–2010), respectively. Most of the projected mean annual snowmelt runoff ratios in third-level basins in different future periods are lower than those in the reference period. Low temperature regions can tolerate more warming, and the snowmelt change in these regions is mainly influenced by precipitation; however, the snowmelt change in warm regions is more sensitive to temperature increases. The spatial variability in snowmelt changes may lead to regional differences in the impact of snowmelt on water supply.

1 Introduction

Snow properties have changed significantly under the ongoing warming of the global climate, and variations in snow cover exert strong feedbacks on the climate system due to its high albedo and low thermal conductivity, as well as the high latent heat of the phase change (Zhang and Ma, 2018; Pulliainen et al., 2020; Vano, 2020; You et al., 2020). Additionally, snow is a critical component of the hydrological system and water cycle, and snowmelt is a major fresh water resource in many regions (Mankin et al., 2015). More than one-sixth of the Earth's population relies on snowmelt for their water supply, and snowmelt-dominated regions contribute roughly one-quarter of the global gross domestic product (Barnett et al., 2005). Climate warming has resulted in smaller snowfall/precipitation ratios and an earlier onset of snowmelt and slower snowmelt rates (Berghuijs et al., 2014; Musselman et al., 2017; Barnhart et al., 2020). This has not only changed seasonal runoff distributions but has also affected the total annual runoff (Bloschl et al., 2019; Jenicek and Ledvinka, 2020). Consequently, determining the amount of snowmelt and its variability under climate change is important for the planning and management of water resources, such as agricultural water management, flood forecasting, reservoir operations and the design of hydraulic structures (Barnhart et al., 2020; Qin et al., 2020).

There are many models for calculating snowmelt, and these can be basically divided into the following two types: physically based snowmelt models and simplified temperature index models (also known as degree–day methods; Skaugen et al., 2018; Li et al., 2019). In theory, physically based snowmelt models can provide more accurate predictions by considering the coupled interaction between energy components in complex snowmelt processes (Li et al., 2019). However, many studies have shown that, due to the mathematical complexities and massive data requirements of physically based models, they do not necessarily perform better than temperature index models (Hock, 2003; Jost et al., 2012; Skaugen et al., 2018; Lopez et al., 2020). Temperature index models are based on the assumption that the temporal variability in incoming solar radiation is adequately represented by the variations in air temperature and that the snowmelt during a time interval is proportional to positive air temperatures (Hock, 2003; Jost et al., 2012; Lopez et al., 2020). Because of the wide availability of air temperature data, their computational simplicity and their generally good model performance, temperature index models are the most common approach for calculating snowmelt around the world (Hock, 2003; Immerzeel et al., 2010; Lopez et al., 2020).

Figure 1 The three main stable snow cover regions and the mean snow depth in China (1951–2009) (a). China's five climatic zones (MPZ – mountain plateau zone; TMZ – temperate monsoon zone; TCZ – temperate continental zone; SMZ – subtropical monsoon zone) and mean annual air temperature (1951–2017) (b). The snow cover classification and mean monthly snow density in China (1999–2008) (c). The third-level basins and mean annual snowfall in China (1961–1979) (d).

China covers a vast area and a variety of climate regions, and its snow-covered regions are widely distributed geographically with evident spatial differences (Tan et al., 2019). Northern Xinjiang, Northeastern China–Inner Mongolia (hereafter referred to as Northeast China) and the Tibetan Plateau are the three main regions with stable snow cover in China (Ke et al., 2016; Fig. 1a). As a typical arid and semi-arid region of Central Asia with a significant lack of freshwater resources, the surface runoff in the Xinjiang Uygur Autonomous Region (hereafter, Xinjiang) is mainly supplied by snowmelt (Chen et al., 2020; Wu et al., 2021). In Northeast China, snow plays an important role as a natural reservoir in winter and a source of water in spring, with snowmelt contributing more than half of the runoff during the main crop sowing season (April and May; Qi et al., 2020). Melting snow is also an important hydrological process on the Tibetan Plateau, which is the source region of many major Asian rivers, and is considered as being the Asian water tower (Immerzeel et al., 2010, 2020). Furthermore, snowmelt is an important water source in other parts of China, especially the arid and semi-arid areas in North China (Zhang et al., 2015; Wu et al., 2021). Spatiotemporal changes in China's snow properties (e.g., snow cover extent, snow cover phenology, snow depth, snow density and snow water equivalent) have drawn wide attention in recent decades (Dai and Chen, 2010; Wang and Li, 2012; Chen et al., 2016; Ke et al., 2016; Tan et al., 2019; Ma et al., 2020; Yang et al., 2020); however, country-wide assessments of snowmelt are lacking. Projected increases in air temperature and associated precipitation regime shifts are expected to have significant consequences for snowmelt and water resources (Ficklin et al., 2016). Although many studies have investigated snowmelt in single or multiple basins (e.g., Zhang et al., 2015, 2020; Chen et al., 2019; Li et al., 2019, 2021), the spatiotemporal variability in snowmelt in China and its response to climate warming remain unclear.

Under this background, we developed a simple monthly temperature index model to calculate the snowmelt in China using a high spatial resolution (0.5^′; approximately 1 km) data set of monthly air temperatures and precipitation. The model considers complex snow processes, such as melting, accumulation and sublimation, and was validated using snowfall, snow depth, snow cover extent and snow water equivalent. The objectives of this study are to (1) quantify the snowmelt across China and in its three main stable snow cover distribution regions, (2) analyze the spatiotemporal variability in snowmelt in China during 1951–2017, (3) elucidate the spatiotemporal variability in snowmelt runoff ratio in the third-level basins in China and (4) assess the impacts of projected future climate change on snowmelt in China.

2 Study region and data sets

2.1 Study region

In general, China can be divided into the following five main climatic zones: the mountain plateau zone (MPZ), the temperate monsoon zone (TMZ), the temperate continental zone (TCZ), the subtropical monsoon zone (SMZ) and the tropical monsoon zone (Yang et al., 2021; Fig. 1b). Because the land area of the tropical monsoon zone is significantly smaller than that of the other climatic zones, and it has extremely low snowfall, it was incorporated into the SMZ for parameter assignment. Additionally, because of the lack of meteorological data, and the fact that it has very little snowfall, Taiwan was not considered in this study. The snow cover types in China can be divided into the following five types: prairie, taiga, tundra, mountain and ephemeral (X. Li et al., 2020; Fig. 1c). The analysis of snowmelt runoff ratio is carried out at the three-level basin scale, and Fig. 1d shows the boundaries of the three-level basins in China (the first-level and second-level basins can be seen in Fig. S1 in the Supplement).

2.2 Data collection

2.2.1 High spatial resolution data set of monthly air temperatures and precipitation

Data with a high spatial resolution of 0.5^′ (approximately 1 km), including the monthly minimum, maximum and mean temperatures (T_min, T_max and T_a, respectively) and precipitation, were obtained from Zenodo in the Network Common Data Form (NetCDF; https://doi.org/10.5281/zenodo.3114194, for precipitation, Peng, 2019a; https://doi.org/10.5281/zenodo.3185722, for air temperature, Peng, 2019b). The data were spatially downscaled from the 0.5^∘ Climatic Research Unit time series data, with the 0.5^′ climatology data set of WorldClim, using delta spatial downscaling. Peng et al. (2019) evaluated this data set using observations collected from 1951 to 2016 by 496 meteorological stations across China, and the average values of the coefficient of determination (R²), mean absolute error (MAE), root mean square error (RMSE) and Nash–Sutcliffe efficiency (NSE), for the monthly mean temperature at all meteorological stations are 0.996, 0.820, 0.969 and 0.981 ^∘C, respectively, and for the monthly precipitation are 0.863, 13.269, 21.941 and 0.808 mm, respectively. Detailed information on the data set was given by Peng et al. (2019). Although the original data set covers the 1901 to 2017 period, we selected 1951 to 2017 as the study period, as only data from after 1951 have been validated by observations from meteorological stations (Table 1).

Table 1 Data used in the snowmelt model.

^* Minimum, maximum and mean temperatures. Note: CMA – China Meteorological Administration; NTPDC – National Tibetan Plateau Data Center; NCDDC – National Cryosphere Desert Data Center.

2.2.2 Data sets for model parameterization

The observational air temperature data used to determine positive degree–day (PDD) parameters were collected from 824 meteorological stations (Fig. 1b) in the 1951–2017 period and were provided by the China Meteorological Administration (http://data.cma.cn/, last access: 27 December 2019). The positive degree–day is defined as the cumulative temperature above 0 ^∘C over a given period (Wake and Marshall, 2017), and in this study, it was the monthly positive accumulated temperature.

The observational snow density data used to determine degree–day factors (DDFs) were collected from 417 meteorological stations in China (Fig. 1c) during the 1999–2008 period and were provided by the China Meteorological Administration.

The threshold temperature parameters for determining the precipitation types at meteorological stations (Fig. S2) were obtained from Han et al. (2010), who derived threshold temperature parameters using air temperature and precipitation data from 485 meteorological stations in China from 1961 to 1979. The threshold temperatures of each calculated cell were interpolated using the parameters from the meteorological stations via the inverse distance weighting (IDW) method.

2.2.3 Data sets for model evaluation

The observational snowfall (snow depth) data used to validate the model were collected from 475 (557) meteorological stations in China (Fig. 1d (1a)) during the 1961–1979 (1951–2009) period and were provided by the China Meteorological Administration (Table 1). Because of the scarce data or the very low snowfall (shallow snow) at some stations, snowfall (snow depth) observations from 457 (264) stations were selected to verify the snowfall (snow depth) from the model.

A long time series of daily snow depth derived from passive microwave remote sensing data (1979–2018) was obtained from the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn, last access: 26 November 2019; Che and Dai, 2015; Table 1). The data set covers the entire land surface of China, with a spatial resolution of 0.25^∘. Detailed information on the data set sources and product processes can be found in Che et al. (2008) and Dai et al. (2012). This data set has been widely utilized in climatic and hydrological research in China (e.g., Liu et al., 2020; Wu et al., 2021; Zhu et al., 2021). The spatial resolution of this data set is significantly different from that of the air temperature and precipitation data used in this study. The snow cover extent measures generated from this data set were used to validate the model.

Additionally, a long time series of a daily snow water equivalent data set derived from passive microwave remote sensing data (1980–2020) was provided by the National Cryosphere Desert Data Center (https://www.ncdc.ac.cn, last access: 9 April 2021; Jiang et al., 2020; Table 1). The snow water equivalent data set was produced from the passive microwave remote sensing data using the mixed-pixel method. The data set covers the entire land surface of China, with a spatial resolution of 25 km. Detailed information on the data set sources and product processes are given by Jiang et al. (2014) and Yang et al. (2019).

2.2.4 Climate projections and downscaling

In total, five CMIP5 models under three different representative concentration pathway scenarios (RCP2.6, RCP4.5 and RCP8.5), namely GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR, MIROC-ESM-CHEM and NorESM1-M, were selected to predict the future snowmelt changes in China. The Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) selected these five CMIP5 models to span the space of global mean temperature change and relative precipitation changes (Warszawski et al., 2014). The aforementioned climate projections have been bias-corrected and downscaled to a grid with a resolution of 0.5^∘ by ISIMIP (Hempel et al., 2013) and widely utilized in climate change research in China (e.g., Li et al., 2016; Yuan et al., 2017; Chang et al., 2020). In this study, we use the delta downscaling method to determine the monthly future meteorological data (2006–2099) based on the high spatial resolution temperature and precipitation data set and the simulations of the five CMIP5 models during the historical period (1951–2005). Detailed information on the delta method was given in Immerzeel et al. (2012) and Zhao et al. (2019). The projected precipitation and temperature changes in China and the three main stable snow cover regions can be seen in Figs. S3 and S4, respectively.

Figure 2 The relationship between the observed monthly accumulated positive air temperature (PDD) and mean air temperature (T_a) (a) and the scatterplot between the calculated and observed monthly PDD (b) at meteorological stations in the mountain plateau zone of China. R² is the coefficient of determination, MAE is the mean absolute error (degrees Celsius), RMSE is the root mean square error (degrees Celsius) and NSE is the Nash–Sutcliffe efficiency.

3 Method

3.1 Snowmelt model

Temperature index models are based on an assumed relationship between snowmelt and air temperature, which is usually expressed in the form of positive temperature sums. In this study, the snowmelt was calculated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{M}_m=min({\text{DDF}}_m\cdot {\text{PDD}}_m\cdot D_m,S_{\text{acc},m})\text{(2)}& {\displaystyle }& {\displaystyle }{S}_{\text{acc},m}=S_{\text{acc},m-\mathrm{1}}+P_{\text{snow},m}-S_{\text{sub},m}-M_m,\end{array}$$

where M, S_acc, P_snow and S_sub are the snowmelt, snow accumulation, snowfall and snow sublimation (millimeters), respectively. DDF is the degree–day factor (millimeters per degree Celsius per day; hereafter $\mathrm{mm}({}^{\circ }\mathrm{C}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$), an empirical factor that relates the rate of snowmelt to air temperature. PDD is the accumulated positive air temperature (degrees Celsius). The subscript m indicates the month, and D is the number of days in the month m.

3.1.1 Snowfall

The determination of the precipitation types is the first crucial step, and the distinction between rainfall and snowfall is based on the assumption that precipitation falls either as rain, as snow or as a mix, depending on two threshold temperature parameters, as follows:

$$\begin{array}{}\text{(3)}& P_{\text{snow}}=\left\{\begin{array}{ll}\mathrm{0}& T_{\mathrm{a}}\ge T_{\text{rain}}\\ P\cdot \frac{T_{\text{rain}}-T_{\mathrm{a}}}{T_{\text{rain}}-T_{\text{snow}}}& T_{\text{snow}}<T_{\mathrm{a}}<T_{\text{rain}}\\ P& T_{\mathrm{a}}\le T_{\text{snow}},\end{array}\right.\end{array}$$

where P_snow and P are the monthly snowfall and precipitation (millimeters), respectively. T_a is the monthly mean temperature (degrees Celsius). T_rain and T_snow are the threshold temperature (degrees Celsius) for rainfall and snowfall, respectively.

3.1.2 DDF

DDF is the key parameter of the temperature index snowmelt model, which depends on the snow density and varies with space and time. We used the following the empirical equation to determine the DDF after Rajkumari et al. (2019):

$$\begin{array}{}\text{(4)}& \text{DDF}=\mathrm{11}({\mathit{\rho }}_{\mathrm{s}}/{\mathit{\rho }}_{\mathrm{w}}),\end{array}$$

where ρ_s and ρ_w are the density of snow and water (grams per cubic centimeter; hereafter g cm⁻³), respectively. The observed snow density is variable and the density of water is constant. DDF is the degree–day factor ($\mathrm{mm}({}^{\circ }\mathrm{C}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$). For taiga regions (Fig. 1c), the equation was corrected as follows:

$$\begin{array}{}\text{(5)}& \text{DDF}=\mathrm{10.4}({\mathit{\rho }}_{\mathrm{s}}/{\mathit{\rho }}_{\mathrm{w}})-\mathrm{0.7}.\end{array}$$

The DDF values at the meteorological stations were calculated for each month and were then interpolated to each calculated cell.

3.1.3 PDD

Based on the daily mean air temperatures from 1960 to 2018 at meteorological stations in China, Liu et al. (2021) proposed a simple method to calculate the monthly PDD according to the monthly mean temperature, as follows:

$$\begin{array}{}\text{(6)}& \text{PDD}=\left\{\begin{array}{ll}\mathrm{0}& T_{\mathrm{a}}\le T_{\mathrm{1}}\\ a\cdot T_{\mathrm{a}}^{\mathrm{2}}+b\cdot T_{\mathrm{a}}+c& T_{\mathrm{1}}<T_{\mathrm{a}}<T_{\mathrm{2}}\\ T_{\mathrm{a}}\cdot D_m& T_{\mathrm{a}}\ge T_{\mathrm{2}},\end{array}\right.\end{array}$$

where T₁ and T₂ are threshold temperatures (degrees Celsius), and a, b and c are empirical parameters.

First, T₁, T₂, a, b and c were calculated according to the measured data from meteorological stations in the four different climate zones (Fig. 1b) and were then interpolated to each calculated cell. Table 2 shows the different threshold temperatures and empirical parameters in the four climate zones, and Fig. 2 shows the relationship between the measured monthly PDD and the mean temperature at meteorological stations and the relationship between the calculated and measured monthly PDD, taking the MPZ as an example. A statistical analysis of the calculated monthly PDD against observed PDD at 824 meteorological stations shows the excellent performance at almost all of the stations (Fig. S5). Statistical analysis shows that Eq. (6) can adequately calculate the monthly PDD (Table 2).

Table 2 The parameters required for the calculation of the monthly accumulated positive air temperature (PDD) and the statistical analysis between the calculated and measured monthly PDD in four different climatic zones of China.

MPZ – mountain plateau zone; TMZ – temperate monsoon zone; TCZ – temperate continental zone; SMZ – subtropical monsoon zone; T₁, T₂, a, b and c – parameters in Eq. (6); R² – coefficient of determination; MAE – mean absolute error (degrees Celsius); RMSE – root mean square error (degrees Celsius); NSE – Nash–Sutcliffe efficiency.

3.1.4 Snow sublimation

As a loss of water from the snowpack to the atmosphere, snow sublimation is difficult to quantify by measurement or modeling (Stigter et al., 2018). The bulk aerodynamic method and aerodynamic profile methods are the common methods to calculate snow sublimation (Svoma, 2016). Some researchers used evapotranspiration equations (e.g., the Penman–Monteith method) to estimate sublimation (Stigter et al., 2018). In general, these methods require accurate meteorological data and are very difficult to scale up from the microscale to the macroscale (Svoma, 2016). Many studies have reported the ratio of snow sublimation to snowfall (e.g., Zhang et al., 2008; Zhu et al., 2014; Sexstone et al., 2018). Therefore, considering the limited data availability and the monthly scale of this study, the following simple equation was used to estimate the monthly snow sublimation:

$$\begin{array}{}\text{(7)}& S_{\text{sub}}=min(k\cdot S_{\text{acc}},\text{PET}),\end{array}$$

where k⋅S_acc and PET are the amount of snow accumulation available for sublimation and the potential sublimation (millimeters), respectively. The empirical parameter k was set according to studies reporting the ratio of snow sublimation to snowfall in China and surrounding areas (Table S1 and Fig. S6), as shown in Table 3.

Table 3 The values of the empirical parameter k used to estimate snow sublimation (Eq. 7) for different climate zones and different snow types.

MPZ – mountain plateau zone; TMZ – temperate monsoon zone; TCZ – temperate continental zone; SMZ – subtropical monsoon zone.

PET was calculated by the Hargreaves–Samani method, as follows (Hargreaves and Samani, 1985):

$$\begin{array}{}\text{(8)}& \text{PET}=\mathrm{0.0023}\cdot R_a(T_{\mathrm{a}}+\mathrm{17.8})(T_{\text{max}}-T_{\text{min}}{)}^{\mathrm{0.5}}/\mathit{\lambda },\end{array}$$

where T_a, T_max and T_min are the monthly average temperature, maximum temperature and minimum temperature (degrees Celsius), respectively; R_a is the extraterrestrial radiation (megajoules per square meter per month; hereafter $\mathrm{MJ}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{month}}^{-\mathrm{1}}$), which is a function of the latitude and the Julian day, and λ is the latent heat of vaporization (megajoules per kilogram; hereafter MJ kg⁻¹).

3.2 Snowmelt runoff ratio

The snowmelt runoff ratio, which represents the contribution of snowmelt to river discharge, is defined as the percentage of runoff derived from snowmelt to the total runoff. Many studies have used indirect indicators to estimate the snowmelt runoff ratio, such as the ratios of total snowfall to total precipitation, total snowfall to total runoff or melt season runoff to total annual runoff (Li et al., 2017). Snowmelt and rainfall are the sources of runoff generation (Vormoor et al., 2016). In this study, the snowmelt runoff ratio was calculated as the ratio of snowmelt to the sum of snowmelt and rainfall, as follows:

$$\begin{array}{}\text{(9)}& {\text{Snow}}_{\mathrm{r}}=M/(M+P_{\text{rain}})\times \mathrm{100},\end{array}$$

where M and P_rain are the snowmelt and rainfall (millimeters), respectively, and Snow_r is the percentage of runoff derived from snowmelt to the total runoff.

3.3 Trend analysis

The non-parametric Mann–Kendall test (Mann, 1945; Kendall, 1975) was used to analyze the trend and significance level of the snowmelt and other variables, as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }Z=\left\{\begin{array}{ll}(S-\mathrm{1})/\sqrt{\text{var}\left(S\right)}& S>\mathrm{0}\\ \mathrm{0}& S=\mathrm{0}\\ (S+\mathrm{1})/\sqrt{\text{var}\left(S\right)}& S<\mathrm{0}\end{array}\right.\text{(11)}& {\displaystyle }& {\displaystyle }S=\sum _{k=\mathrm{1}}^{n-\mathrm{1}}\sum _{j=k+\mathrm{1}}^n\text{sgn}(x_j-x_k)\text{(12)}& {\displaystyle }& {\displaystyle }\text{sgn}(x_j-x_k)=\left\{\begin{array}{ll}\mathrm{1}& (x_j-x_k)>\mathrm{0}\\ \mathrm{0}& (x_j-x_k)=\mathrm{0}\\ -\mathrm{1}& (x_j-x_k)<\mathrm{0}\end{array}\right.\text{(13)}& {\displaystyle }& {\displaystyle }\mathrm{\text{var}}\left(S\right)={\displaystyle \frac{n\left(n-\mathrm{1}\right)\left(\mathrm{2}n+\mathrm{5}\right)-{\sum }_{i=\mathrm{1}}^m{t}_i(t_i-\mathrm{1})(\mathrm{2}{t}_i+\mathrm{5})}{\mathrm{18}}},\end{array}$$

where n is the number of the data set, x_j and x_k are the data values in time series j and k, m is the number of tied groups, and t_i denotes the number of ties of extent i. A tied group is a set of sample data having the same value. Positive and negative values of Z indicate increasing and decreasing trends, respectively. The testing of the trends is done at the specific α significance level. If $\left|Z\right|>Z_{\mathrm{1}-\mathit{\alpha }/\mathrm{2}}$, then the trend is statistically significant; otherwise, the trend is not statistically significant. In this study, significance levels of α=0.05 (95 % confidence level) were applied.

Additionally, Sen's (1968) slope method was used to analyze the slope of the variation, as follows:

$$\begin{array}{}\text{(14)}& \mathit{\beta }=\text{median}\left({\displaystyle \frac{x_j-x_k}{j-k}}\right)j>k,\end{array}$$

where the β sign reflects whether a trend is negative or positive, while its value indicates the steepness of the trend.

3.4 Evaluation criteria

Statistical indices were used for quantitative analysis of the snowmelt model performance, and a series of statistical criteria are used in this study as follows:

$$\begin{array}{}\text{(15)}& R^{\mathrm{2}}={\left[\text{cov}(X_{\mathrm{s}},X_{\mathrm{o}})/\mathit{\sigma }{X}_{\mathrm{s}}\mathit{\sigma }{X}_{\mathrm{o}}\right]}^{\mathrm{2}}\end{array}$$

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }\text{MAE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left|X_{\mathrm{s}i}-X_{\mathrm{o}i}\right|\text{(17)}& {\displaystyle }& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(X_{\mathrm{s}i}-X_{\mathrm{o}i}\right)}^{\mathrm{2}}}{n}}}\text{(18)}& {\displaystyle }& {\displaystyle }\text{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(X_{\mathrm{s}i}-X_{\mathrm{o}i}\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(X_{\mathrm{o}i}-\stackrel{\mathrm{‾}}{X_{\mathrm{o}}}\right)}^{\mathrm{2}}}},\end{array}$$

where R², MAE, RMSE and NSE are the coefficient of determination, mean absolute error, root mean square error and Nash–Sutcliffe efficiency, respectively. X_si and X_oi represent simulated and observed data at time i, respectively. n is the number of data points, and cov and σ are the covariance and standard deviation, respectively.

4 Results

4.1 Model performance

4.1.1 Snowfall validation using observational data from meteorological stations

Figure 3 shows the statistical criteria between the calculated and observed snowfall at meteorological stations in China from 1961 to 1979. The MAE and RMSE vary greatly between stations due to the huge difference in the amount of snowfall at each station. The ratios of MAE and RMSE to monthly mean snowfall ($\text{MAE}/\text{Mean}$ and $\text{RMSE}/\text{Mean}$, respectively) were selected to analyze the model performance. Of the 457 stations, 315 (68.9 %) have R²>0.4, 263 (57.5 %) have $\text{MAE}/\text{Mean}<\mathrm{1}$, 274 (60.0 %) have $\text{RMSE}/\text{Mean}<\mathrm{1}$, and 353 (77.2 %) have NSE>0. Among China's three main stable snow cover regions, the most accurate snowfall simulation is obtained for Northeast China, followed by North Xinjiang and then the Tibetan Plateau.

Figure 3 Statistical criteria of the calculated snowfall against observed snowfall at 457 meteorological stations in China. Panel (a) shows the R² (coefficient of determination). Panel (b) shows $\text{MAE}/\text{Mean}$, where MAE is the mean absolute error, and Mean is the monthly mean snowfall. Panel (c) shows $\text{RMSE}/\text{Mean}$, where RMSE is the root mean square error. Panel (d) shows the NSE (Nash–Sutcliffe efficiency).

4.1.2 Snow depth validation using observational data from meteorological stations

The monthly snow depth in each grid was calculated from the snow accumulation and snow density. As the snow depth was observed at meteorological stations daily, the snow depth observation on the last day of each month was selected to verify the model's snow depth output. Figure 4 shows the statistical analyses of the calculated and observed snow depths in China from 1951 to 2009. Similar to snowfall validation, besides R² and NSE, the ratios of MAE and RMSE to the mean snow depth ($\text{MAE}/\text{Mean}$ and $\text{RMSE}/\text{Mean}$, respectively) were selected to analyze the model performance. Of the 264 stations, 108 (40.9 %) have R²>0.2, 221 (83.7 %) have $\text{MAE}/\text{Mean}<\mathrm{2}$, 52 (19.7 %) have $\text{RMSE}/\text{Mean}<\mathrm{2}$, and 105 (39.8 %) have NSE>0. Better snow depth simulations are obtained in regions with larger snow depths, such as Northeast China and North Xinjiang.

Figure 4 Statistical criteria of the calculated snow depth against observed snow depth at 264 meteorological stations in China. Panel (a) shows the R² (coefficient of determination). Panel (b) shows $\text{MAE}/\text{Mean}$, where MAE is the mean absolute error, and Mean is the monthly mean snow depth. Panel (c) shows $\text{RMSE}/\text{Mean}$, where RMSE is the root mean square error. Panel (d) shows the NSE (Nash–Sutcliffe efficiency).

4.1.3 Snow cover extent validation using the data set derived from passive microwave remote sensing

As the extent of the snow cover derived from remote sensing has a daily timescale and the snowmelt model has a monthly timescale, the remote sensing data on the last day of each month were selected to verify the model's snow cover extent output. The mean monthly snow cover extent in China output by the model and the information derived from passive microwave remote sensing from 1979 to 2017 are 1.74×10⁶ and 1.62×10⁶ km², respectively. The R², MAE, RMSE and NSE between the snow cover extent output by the model against the data set derived from remote sensing in China from 1979 to 2017 are 0.93, 0.45×10⁶, 0.64×10⁶ km² and 0.89, respectively. Among the three main stable snow cover regions, the R², MAE, RMSE and NSE are, respectively, 0.81, 0.06×10⁶, 0.09×10⁶ km² and 0.76, for Northern Xinjiang, 0.93, 0.13×10⁶, 0.21×10⁶ km² and 0.87 for Northeast China, and 0.90, 0.32×10⁶, 0.40×10⁶ km² and 0.81 for the Tibetan Plateau (Fig. 5).

Figure 5 Scatterplots of the snow cover extent output by the model and information derived from remote sensing in China and its three main stable snow cover regions. The red and blue solid lines are the linear fit and the 1:1 line, respectively.

4.1.4 Snow water equivalent validation using the data set derived from passive microwave remote sensing

The snow water equivalent derived from passive microwave remote sensing on the last day of each month was used to verify the snow accumulation output by the model. As the spatial resolution of the snow water equivalent data set is 25 km, and the grid scale of the model is about 1 km, to facilitate comparison, they were uniformly converted into water equivalent units of cubic meters for different regions. The mean monthly snow accumulation in China output by the model from 1980 to 2017 is 2.55×10¹⁰ m³, while the mean snow water equivalent derived from passive microwave remote on the last day of each month is 1.47×10¹⁰ m³. The R², MAE, RMSE and NSE between the snow accumulation output by the model against snow water equivalent data set in China from 1980 to 2017 are 0.62, 1.27×10¹⁰, 1.27×10¹⁰ m³ and 0.80, respectively. The R², MAE, RMSE and NSE for the three main stable snow cover regions are shown in Fig. 6. A large number of glaciers are distributed in the Tibetan Plateau, many of which may not be recorded in the remote-sensing-based snow water equivalent data set because its spatial resolution of about 25 km is larger than the width of most glaciers. This can explain the fact that a snow water equivalent of 0 is observed in some months in this region. However, the spatial scale of the model is about 1 km, meaning that the snow accumulation in the glacier areas is always detected. This can explain the observation that, in the Tibetan Plateau, the modeled snow water equivalent is higher than the remote-sensing-based values. When the snow accumulation in the glacier areas was not considered, the model performance improved in both the Tibetan Plateau and the whole of China (MAE and RMSE decreased and NSE increased; see Fig. S7).

Figure 6 Scatterplots of the snow water equivalent output by the model and derived from remote sensing in China and its three main stable snow cover regions. The red and blue solid lines are the linear fit and the 1:1 line, respectively.

4.2 Spatial and temporal variability in snowmelt

4.2.1 Spatial distribution

The mean annual snowmelt in China from 1951 to 2017 is about 2.41×10¹¹ m³ yr⁻¹. The three main stable snow cover regions account for about 72.62 % of the total snowmelt in China, with the mean annual snowmelt in Northern Xinjiang, Northeast China and the Tibetan Plateau being about 0.18×10¹¹, 0.42×10¹¹ and 1.15×10¹¹ m³ yr⁻¹, respectively. Southeast China (mainly the Huaihe River basin and the lower reaches of the Yangtze River) also has high snowmelt due to heavy snowfall (Fig. 7a). The areas with the lowest snowmelt are mainly distributed in the arid region of Northwest China, because of this region's low precipitation, and in the humid region of South China, because of this region's high air temperature and low snowfall.

Figure 7 Spatial distribution of mean annual snowmelt (a), the month of maximum snowmelt (b), trends in annual snowmelt based on the Mann–Kendall method (c) and Sen's slope of the annual snowmelt in China during the 1951–2017 period.

Figure 8 Spatial distribution of the mean monthly snowmelt (millimeters) for 12 months in China during the 1951–2017 period.

The spatial pattern of snowmelt differs considerably during the year (Fig. 8). In the three stable snow cover regions, in winter, the air temperature is low, and there is little snowmelt. The winter snowmelt mainly occurs in South China. In spring, North Xinjiang and Northeast China are the main regions of snowmelt in China. With further warming, the Tibetan Plateau becomes the main region of snowmelt until May. In summer, there is no snowfall in most of China, and snowmelt occurs only in the high mountains of the Tibetan Plateau and the Tianshan mountains.

Figure 9 Interannual variability in the mean annual snowmelt in China and its three main stable snow cover regions from 1951 to 2017.

The month of maximum snowmelt differs between regions (Fig. 7b). In Northern Xinjiang and Northeast China, due to the warming in spring, the maximum monthly snowmelt generally occurs in March, April or May. Meanwhile, in North and Southeast China, which are ephemeral snow regions where the snow melts quickly after falling, the maximum snowmelt generally occurs in the month with the largest snowfall (which occurs in winter). Because of the Tibetan Plateau's complex terrain and varied climate, the months of maximum snowmelt vary greatly across this region. In the southeastern part of the Tibetan Plateau, the month of maximum snowmelt occurs during the winter snowfall period because of the warm and humid climate. In the area between the Himalayas and the Gangdise Mountains, the Qaidam Basin and other colder regions (Fig. S8), the month of maximum snowmelt generally occurs in spring. Meanwhile, in the Qiangtang plateau, which has a drier climate, most snowfall occurs in summer, and the month of maximum snowmelt also occurs in this season. There are many high-elevation mountains in the Tibetan Plateau, with air temperatures above 0 ^∘C for only a few days in the summer, and the month of maximum snowmelt in those high elevations areas also occurs in summer.

4.2.2 Temporal variations

From 1951 to 2017, the annual snowmelt increased significantly in some areas of Northern Xinjiang, although it did not increase significantly in Northern Xinjiang as a whole. In Northeast China, although the central and eastern regions showed a significant increasing trend, the southwestern regions showed a decreasing trend, leading to a slight decreasing trend (although not a significant trend) in the whole of Northeast China. Annual snowmelt increased significantly in most parts of the Tibetan Plateau, leading to a significant increasing trend in the Tibetan Plateau as a whole. In Southeast China, the annual snowfall decreased due to climate warming, and snowmelt also decreased significantly in this region, leading to a decreasing trend (although not a significant trend) of annual snowmelt in China (Figs. 7c and d and 9). From 1951 to 2017, the linear trends of annual snowmelt in China, Northern Xinjiang, Northeast China and the Tibetan Plateau are $-\mathrm{2.7}\times {\mathrm{10}}^{\mathrm{9}}$, 0.2×10⁹, $-\mathrm{0.3}\times {\mathrm{10}}^{\mathrm{9}}$ and 1.4×10⁹ m³ per decade, respectively.

The temporal trend of snowmelt in each month from 1951 to 2017 varies across China (Figs. S9 and S10). In Northern Xinjiang and Northeast China, snowmelt increased significantly in March but decreased significantly in April. This might imply that the warming of these two regions led to the earlier onset of snowmelt. Regarding seasonal changes between 1951 and 2017, at the Tibetan Plateau, spring snowmelt increased significantly, whereas summer snowmelt decreased significantly. This may be due to the warming in spring and the reduction of snowfall in summer. In Southeast China, snowmelt experienced a decreasing trend in almost all months because of the reduction in snowfall due to climate warming.

4.3 Spatial and temporal variability in the snowmelt runoff ratio

4.3.1 Spatial distribution

Of the 210 third-level basins in China, only nine small basins near the tropical monsoon region have no snowmelt runoff from 1951 to 2017 (as snowmelt was not calculated in Taiwan, it was assumed that there was no snowmelt there). In West China, which contains two main stable snow cover regions – Northern Xinjiang and the Tibetan Plateau – the mean annual snowmelt runoff ratios are more than 10 % in all basins, except those in the Taklimakan Desert. In the basins in North and Northeast China, the snowmelt runoff ratios are generally more than 5 %, whereas, due to heavy rainfall, the snowmelt runoff ratios in basins in South China are generally less than 2 % (Fig. 10a).

Figure 10 Spatial distribution of the mean annual snowmelt runoff ratio (percent) (a), and Sen's slope of the annual snowmelt runoff ratio (percent per year; hereafter % yr⁻¹) (b) in the third-level basins in China during the 1951–2017 period.

Figure 11 Spatial distribution of the monthly mean snowmelt runoff ratio (percent) for 12 months in the third-level basins in China during the 1951–2017 period.

The monthly snowmelt runoff ratio also shows large spatial differences (Fig. 11). In the cold months of November, December, January and February, the snowmelt runoff ratios in the basins of Central and North China are over 30 % due to the precipitation being dominated by snowfall because of the extremely low temperatures in the cold months, and there is no snowmelt or rainfall in some basins of the three stable snow cover regions. In March, although the rainfall increased, the snowmelt from the snow that accumulated in winter also increased, and the snowmelt runoff ratios are still relatively high in most basins of North China. In April, the snowmelt runoff ratios begin to decline in most basins, and basins with snowmelt runoff ratios greater than 30 % are mainly located in the three stable snow cover regions. In May and June, the snowmelt runoff ratios further decrease, dropping to zero in most of South, East, and Northeast China; however, they are still more than 30 % in some basins in Northern Xinjiang and the Tibetan Plateau. In July and August, there is no snowmelt runoff in any basins, except for some in Xinjiang and the Tibetan Plateau, and almost no basins have a snowmelt runoff ratio of more than 30 %. In September and October, snowfall and snowmelt gradually increase, leading to a gradual increase in the snowmelt runoff ratio in many basins.

4.3.2 Temporal variations

From 1951 to 2017, the Sen's slope in third-level basins shows that the annual snowmelt runoff ratio decreased in most basins (Fig. 10b). Additionally, the Mann–Kendall test shows that the basins with a significant decreasing trend of the snowmelt runoff ratio are mainly distributed in central Inner Mongolia, the southern slope of the Tianshan Mountains and South China (Fig. S11). The basins with an increased snowmelt runoff ratio are mainly distributed in the southeastern part of the Tibetan Plateau, the Heihe River basin in the Qilian Mountains, the Gurbantünggüt Desert and Wulungu River in Northern Xinjiang, and the Songhua River basin in Northeast China. Among these basins, the Mann–Kendall test shows that only four basins had significant increasing trends, i.e., the source region of the Yellow River and three subbasins of the Songhua River.

The temporal trend in the snowmelt runoff ratio from 1951 to 2017 shows spatial variations in China in every month (Figs. S12 and S13). In the third-level basins in Central and East China, the snowmelt runoff ratio decreased in almost every month. In December, January and February, the snowmelt runoff ratio decreased significantly in the subbasins of the Huaihe River and the middle and lower reaches of the Yangtze River. In March, the snowmelt runoff ratio decreased mainly in the middle reaches of the Yellow River and in Northwest China. In April and May, the snowmelt runoff ratios decreased significantly in Northern Xinjiang and Northeast China, and in June, July and August, they decreased significantly in the Tibetan Plateau, Tianshan mountains and Altai Mountains. In September, the snowmelt runoff ratios decreased in only a few basins, mainly in the Tianshan mountains and at the edge of the Tibetan Plateau. In October, the snowmelt runoff ratios decreased significantly, mainly in Northern Xinjiang and central Inner Mongolia in October, and in November, they decreased significantly in Southern Xinjiang and the middle reaches of the Yellow River. In some months, the snowmelt runoff ratio increased in a few basins, mainly in the Tibetan Plateau and nearby areas. For example, in May, the snowmelt runoff ratios showed a large increase in the upper reaches of four rivers, namely the Yangtze River, the Nu River, the Lancang River and the Yarlung Zangbo River; however, the Mann–Kendall test shows that these monthly increasing trends in third-level basins are rarely significant (Fig. S13).

4.4 Future changes in snowmelt under different climate scenarios

4.4.1 Snowmelt

The snowmelt in Northern Xinjiang and Northeast China show a significant decreasing trend from 2006 to 2099 under RCP8.5 but shows no significant changes under RCP2.6 or RCP4.5. The snowmelt in the Tibetan Plateau and the whole of China show significant decreasing trends from 2006 to 2099 under all three RCPs, with the most drastic decline being under RCP8.5 (Fig. 12). The changes in snowmelt in China and its three main stable snow cover regions in the near-future (2011–2040), mid-future (2041–2070), far future (2071–2099) periods under the three RCPs are shown in Table 4. In Northern Xinjiang, the total projected snowmelt in these three future periods under the three RCPs are not very different than the snowmelt in the reference period (1981–2010). The models project that snowmelt will increase in low-elevation arid areas and decrease in the higher-elevation Tianshan and Altai Mountains (Fig. 13). In Northeast China, the total projected snowmelt in different future periods under the three RCPs is lower than in the reference period. The models project that snowmelt will increase in the Greater Khingan Range and the Songliao Plain and decrease in the Lesser Khingan and Changbai mountains. In most areas of the Tibetan Plateau, the models project a large decrease in snowmelt (Fig. 13). Under RCP2.6, RCP4.5 and RCP8.5, the snowmelt in the Tibetan Plateau in the near future (mid-future; far future) is projected to decrease by 16.7 % (20.8 %; 19.2 %), 18.8 % (24.8 %; 26.1 %) and 18.3 % (31.1 %; 44.2 %), respectively, compared to the reference period. Southeast China is another area where snowmelt is projected to decrease to a large degree (Fig. 13). The models project that the total snowmelt in China will decrease in different future periods. Under RCP2.6, RCP4.5 and RCP8.5, the projected decrease in snowmelt in China in the near future (mid-future; far future) is 10.4 % (15.8 %; 13.9 %), 12.0 % (17.9 %; 21.1 %) and 11.7 % (24.8 %; 36.5 %), respectively, compared to the reference period.

Figure 12 Projected future changes in snowmelt in China and its three main stable snow cover regions. The historical period is 1951–2017, and RCP is the representative concentration pathway.

Table 4 Changes in mean annual air temperature (T_a; degrees Celsius), precipitation (P; millimeters) and snowmelt (M; ×10¹¹ m³) in China and its three main stable snow cover regions during the reference period (1981–2010), near future (2011–2040), mid-future (2041–2070) and far future (2071–2099) under the three representative concentration pathways (RCPs).

Figure 13 Spatial distributions of differences (millimeters) between the projected mean annual snowmelt and the reference period (1981–2010) for the near future (2011–2040), mid-future (2041–2070) and far future (2071–2099) under RCP2.6, RCP4.5 and RCP8.5 in China.

4.4.2 Snowmelt runoff ratio

Under the three RCPs, the projected mean annual snowmelt runoff ratios in the third-level basins in different future periods are mostly smaller than those in the reference period, except for a few basins in Xinjiang and North China (Fig. 14). In general, the largest decreases in snowmelt runoff ratio in the basins are projected to occur by in far future, followed by the mid-future and near future. The largest decreases are projected under RCP8.5, followed by RCP4.5 and RCP2.6. Among the three main stable snow cover regions, the snowmelt runoff ratios are projected to decrease the most in basins in the Tibetan Plateau, followed by basins in Northern Xinjiang and Northeast China. Under RCP8.5, the projected mean annual snowmelt runoff ratios in the far future are lower than those in the reference period in all basins, except the three basins near the Taklimakan Desert and one basin in central Inner Mongolia. Under RCP8.5, relative to the reference period, the snowmelt runoff ratios in the Tibetan Plateau and Tianshan mountains are projected to decrease by more than 5 % in most basins and by more than 10 % in a few basins in the far future (Fig. 14i).

Figure 14 Spatial distributions of differences (percent) between the projected mean annual snowmelt runoff ratio and the reference period (1981–2010) for the near future (2011–2040), mid-future (2041–2070) and far future (2071–2099) under RCP2.6, RCP4.5 and RCP8.5 in China.

5 Discussion

5.1 Model evaluation

Snowmelt is difficult to measure directly, and therefore, the model outputs of snowfall, snow depth, snow cover extent and snow water equivalent were selected to verify the model performance. Although the model is solely driven by air temperature and precipitation, the model outputs show acceptable performance compared with the results from other studies. Han et al. (2021) used rainfall and snowfall temperature thresholds to identify precipitation types in the Lancang River basin in Southwest China and obtained R² values between the simulated snowfall and the snowfall observed at three meteorological stations of 0.42, 0.34 and 0.61, respectively; in this study, the R² values at the same stations are 0.72, 0.46 and 0.39, respectively. Z. Li et al. (2020) used temperature thresholds to calculate the snowfall in the Tianshan mountains of Central Asia and obtained a mean R² value between the simulated snowfall and the snowfall observed at 27 meteorological stations of 0.61; in this study, the mean R² value at 50 stations in Xinjiang is 0.39. Zhong et al. (2018) discriminated the precipitation phase based on temperature thresholds in the Songhua River basin, Northeast China, and obtained R² values between simulated and observed snowfall of less than 0.3 for most meteorological stations; in this study, the R² is larger than 0.3 for most stations.

Snow depth is the worst-performing model output, and this is mainly for the following reasons. (1) The model was performed on a monthly scale, whereas the output data were compared with the observed snow depth measured on the last day of each month, which increased the error. (2) The output data are at the grid scale, whereas the observed snow depth is at the site scale, where the snow properties are not always representative of snow at grid scale (Sexstone et al., 2020). (3) Snow depth itself is difficult to simulate, and even the snow depth retrieved by remote sensing has been shown to have high uncertainty. Orsolini et al. (2019) found that the mean annual snow depth at 33 meteorological stations in the Tibetan Plateau, based on multiple global reanalysis products, was 1.38 to 11.71 cm, with a mean value of 7.88 cm, while the observed snow depth was 0.23 cm. Furthermore, Bin et al. (2013) evaluated snow depth obtained from five algorithms using Advanced Microwave Scanning Radiometer for EOS (AMSR-E) passive microwave against ground observations from meteorological stations across China and found that the RMSE varied from 6.85 to 16.79 cm in Xinjiang region and from 6.21 to 18.05 cm in Northeast China. In this study, the RMSE varies from 0.56 to 13.32 cm in Xinjiang and 0.54 to 9.00 cm in Northeast China. Additionally, many studies have shown that the retrieved snow depth is more accurate in regions with larger snow depth (e.g., Zhou et al., 2017; Wang and Zheng, 2020), and in this study, the performances in the regions with larger snow depth such as Northern Xinjiang and Northeast China are also much better than those in other regions.

Among the simulated snow properties, the snow cover extent shows the best performance, with R² and NSE values above 0.80 (Fig. 5). The performance of the snow water equivalent is acceptable in northern Xinjiang and Northeast China but is poor in the Tibetan Plateau (Fig. 6). There are several possible reasons for this difference. (1) The accuracy of the driving data of precipitation and temperature in the Tibetan Plateau is lower than that in other regions (Peng et al., 2019). (2) Due to the sparse distribution of meteorological stations in the Tibetan Plateau and the fact that most of these are located at low elevations, the reliability of the model parameters might be worse in this region than in other regions. (3) The snow water equivalent data used for verification have large uncertainties in the Tibetan Plateau because of its high elevation and complex terrain and climatic conditions (You et al., 2020).

In summary, the verification of snowfall, snow depth, snow cover extent and snow water equivalent suggests that the model is reliable for calculating snowmelt in China.

5.2 Influence of changes in temperature and precipitation on snowmelt

Snowmelt is sensitive to both temperature and precipitation, and the relationship between snowmelt and warming is more complex than monotonic declines (Mankin et al., 2015). In China, climate warming has led to temperature increase, and snowmelt is increasing significantly in some regions (Fig. 7c). The grids with significant changes in snowmelt from 1951 to 2017 were selected to analyze the changes in air temperature, precipitation and snowfall (Fig. 15). In the regions where snowmelt changed significantly, the temperature generally increased (although not all significantly). Precipitation decreased (increased) significantly in many regions where snowmelt increased (decreased) significantly. The regions with significant changes in snowmelt are relatively consistent with the regions with snowfall changes. There are no grids where snowfall increased significantly and snowmelt decreased significantly. In a few grid cells, snowfall decreased significantly and the snowmelt increased significantly. These grids are too few and can be ignored (Table S2). Compared with the total annual precipitation, the change in snowfall has more influence on snowmelt. Tan et al. (2019) analyzed the relationship between snow cover days and precipitation in China and found that snow cover days were highly correlated with winter precipitation but were not correlated with spring precipitation. In global climate models, there is more uncertainty about precipitation than temperature (Woldemeskel et al., 2016). Many assessments of precipitation from climate models were performed at the annual scale (e.g., Woldemeskel et al., 2016; Ahmed et al., 2019; Yue et al., 2019). Assessing precipitation from climate models using multiple timescales or seasons may allow a more accurate prediction of snowfall and snowmelt changes.

Figure 15 Trends of the annual snowmelt, precipitation, snowfall and air temperature based on the Mann–Kendall method in China during the 1951–2017 period. The meanings of the numbers in the legend are shown in Table S2.

As snow cover formation and snowmelt are closely related to the temperature threshold of 0 ^∘C, temperature change near 0 ^∘C is more likely to trigger a drastic change in snowmelt. This may partly explain the spatial distribution of snowmelt changes in China. In Northern Xinjiang and Northeast China, the main stable snow cover regions in China, the changes in snowmelt from 1951 to 2017 were not significant, and the future projections suggest that the changes in snowmelt in these regions will not be drastic. In these two regions, the snow season temperatures are well below freezing point (Qin et al., 2006), and the snowmelt trend depends on the snowfall trend rather than the temperature trend. There is no reason for snowmelt to change significantly when precipitation does not change significantly and temperature is still to remain well below freezing in those two regions. Therefore, in Northern Xinjiang and Northeast China, snowmelt is unlikely to change significantly when precipitation does not change significantly and the temperature remains well below freezing. In Southeast China, the snowmelt decreased significantly from 1951 to 2017, and the model projects that snowmelt will decrease further under RCP2.6, RCP4.5 and RCP8.5. The climate in this region is warm, with average monthly temperatures near 0 ^∘C in winter (Zeng et al., 2016). Precipitation types are sensitive to temperature increase, and under climate change, less precipitation falls as snow and snowmelt decreases. With the high elevation and low temperature, the Tibetan Plateau experienced the most complex changes in snowmelt from 1951 to 2017. Temperature significantly increased in the Tibetan Plateau, and the annual precipitation increased in most areas of the region (Kuang and Jiao, 2016). This study shows that the snowmelt in the Tibetan Plateau increased significantly from 1951 to 2017. High elevations tend to be colder and therefore can tolerate more warming (Livneh and Badger, 2020). The increased precipitation in the Tibetan Plateau may have offset the decrease in snowfall caused by the temperature increase, and precipitation may be a more influential control on snowmelt in that region than temperature. The results of the CMIP5 models show that the snowmelt in the Tibetan Plateau may decrease in the future. In this region, further warming could cause the temperature threshold for snow formation to be reached, resulting in less snowfall and snowmelt at medium and low elevations, and the temperature may become a more influential control on snowmelt.

5.3 Impact of snowmelt change on regional water supply

China's three main stable snow cover regions (Northern Xinjiang, Northeast China and the Tibetan Plateau) are located in three different climate areas (TCZ, TMZ and MPZ, respectively). Besides these three regions, the part of Southeast China (mainly the Huaihe River basin and the lower reaches of the Yangtze River) that is in the SMZ also has relatively high snowmelt (Fig. 7a). In the following, we examine these four regions to shed light on the potential regional water supply problems that may be associated with changes in snowmelt under climate change.

5.3.1 Northern Xinjiang

Northern Xinjiang is located in the TCZ, and snowmelt is an important source of fresh water in this region (Chen et al., 2020; Wu et al., 2021). This study shows that between 1951 and 2017 the mean annual snowmelt runoff ratios in the third-level basins in Northern Xinjiang are above 10 % and exceeded 30 % in some basins. Although the annual snowmelt showed a significant increasing trend from 1951 to 2017, the increase was mainly in November. Because of the earlier onset of snowmelt under climate warming, the monthly snowmelt increased in March and decreased significantly in April. The snowmelt runoff ratios in April are generally more than 30 % and also show a significant decreasing trend. March and April are times of plant and crop cultivation, and the water requirements for agricultural irrigation therefore increase sharply in these months. The shift of the snowmelt time may reduce the water that is available for agriculture because of the mismatch between snowmelt and crop growing (Notarnicola, 2020) and, therefore, introduce agricultural risks to the Northern Xinjiang.

5.3.2 Northeast China

In Northeast China, the runoff in March and April is mainly generated from snowmelt. From 1951 to 2017, snowmelt in March increased while that in April decreased significantly. The snowmelt runoff ratios in the basins in Northeast China are generally about 10 %, and the importance of snowmelt as a water resource in this region may be less than in Northern Xinjiang. However, because of the cold and long winter, snow plays a role as a natural reservoir to store water in winter, and melting snow releases water in spring, which can influence agriculture by affecting soil moisture (Qi et al., 2020). Snowmelt is vital for the seasonal water supply in Northeast China, which is an important region for agricultural production in China. Because of the changes in snowmelt documented in this study, changes in water resource allocation in spring should be considered to meet the demands of spring crop planting.

5.3.3 Tibetan Plateau

The Tibetan Plateau is known as the water tower of Asia, since many of the continent's major rivers originate there (e.g., the Yangtze, Yellow, Indus, Mekong, Brahmaputra, Salween and Ganges rivers). The western and northern parts of the Tibetan Plateau are the source of many arid inland rivers such as the Hotan and Heihe rivers. Additionally, the Tibetan Plateau contains arid inland basins (e.g., the Qiangtang inland basin). These waters from the Tibetan Plateau have been sustaining life, agricultural, and industrial water usage for nearly 40 % of the world's population (Xu et al., 2008). Due to its high elevation and low temperature, the mean annual snowmelt runoff ratios in the basins of the Tibetan Plateau are generally greater than 10 %. Snowmelt contributes to runoff every month, and in some months, there is no rainfall to generate runoff. The snowmelt in the Tibetan Plateau showed an increasing trend from 1951 to 2017, and as the temperature is projected to rise further in the future, it is likely that both the snowmelt and snowmelt runoff ratio will decrease in the near future. Other studies have similarly concluded that the snowmelt and snowmelt runoff ratio in the Tibetan Plateau will decline in the future (Qin et al., 2020; Kraaijenbrink et al., 2021). The Tibetan Plateau and its surrounding downstream areas are projected to experience declines in the share of water from snowmelt and will thus require increases in alternative water supplies.

5.3.4 Southeast China

Southeast China is located in the SMZ, where precipitation is relatively high. The snowfall in winter is high, resulting in a large amount of snowmelt in this region. However, due to this region's high precipitation, snowmelt contributes relatively little to its water resources, with the snowmelt runoff ratios being less than 2 % in most basins. Although climate change has caused a significant decrease in snowmelt in this region, and this trend is projected to continue in the future, the impact on the water supply is likely to be much smaller than in other regions of China. The increased frequency and intensity of flooding that is projected to be caused by climate change is a major concern for this region (Qin and Lu, 2014).

5.4 Uncertainties and limitations

The snowmelt simulation reported in this study was largely driven by the climatology data set of WorldClim, using delta spatial downscaling, and any uncertainties in the data set were likely to propagate to the snowmelt simulation. Many studies showed that, although the WorldClim data were adequate for air temperatures, there were large errors in its precipitation data, especially in mountainous areas with complex topography and high elevations (Beck et al., 2020; Bobrowski et al., 2021). Although the downscaled data set generated by the bilinear interpolation method containing detailed topographic information, and the effects of distance to the nearest coast and satellite-derived covariates, is of high quality, there is still a gap between the downscaled data and the observed data, especially in the mountainous areas (Peng et al., 2019; Ding and Peng, 2020). Uncertainties in the CMIP5 models (Woldemeskel et al., 2016) and the downscaling method (Lima et al., 2021) may introduce uncertainties in the assessments of future changes in snowmelt.

As the results of this study are based on simple temperature index model simulations, the uncertainty arising from the model parameterization needs to be addressed. There are three important parameters in the model, namely the threshold temperatures for the separation of rainfall and snowfall (T_rain and T_snow), the degree–day factors (DDFs) and the ratio of snow sublimation to snow accumulation (k). Although T_rain and T_snow were derived from observations at meteorological stations, errors and uncertainties may have been increased when the values were interpolated into the grid. The DDF was estimated using an empirical method that depends on the snow density observed at meteorological stations. Due to the uneven distribution of meteorological stations, especially few and sparse meteorological stations in the Tibetan Plateau, and the uncertainties introduced by the spatial interpolation of data from meteorological stations (Zhang and Ma, 2018), the errors of these parameters may have led to the uncertainty in the snowmelt calculation. Additionally, snow sublimation is very difficult quantify by measurement or modeling (Stigter et al., 2018). Due to the limited data availability and the complexity of snow sublimation, this study used the ratio of snow sublimation to snow accumulation to simplify the calculation of the sublimation. However, this ratio can vary considerably, both spatially and temporally, having been estimated to vary between 0.1 % and 90 % around the world (Stigter et al., 2018). The simplified calculation of snow sublimation and the failure to consider the monthly variability of this parameter may have increased the uncertainty of the snowmelt calculation.

Moreover, there are many glaciers in West China, and there are differences in the processes and parameters of glacier melting and snow melting (Terink et al., 2015; Armstrong et al., 2018). In this study, snowmelt was not clearly distinguished in the glacier area, which is another source of uncertainty in the snowmelt calculation.

In addition to snowmelt, the snowmelt runoff ratio is another important research object in this study. We estimated this ratio using $\text{snowmelt}/(\text{snowmelt}+\text{rainfall})$. However, this method may underestimate the ratio because snowmelt is more effective at generating catchment runoff than rainfall (Li et al., 2017; Jenicek and Ledvinka, 2020).

6 Conclusions

A simple temperature index model was used to calculate a time series of monthly snowmelt at a resolution of 0.5^′ (approximately 1 km) in China for the 1951–2017 period. The mean annual snowmelt in China is 2.41×10¹¹ m³ yr⁻¹, and the mean annual snowmelt in Northern Xinjiang, Northeast China and the Tibetan Plateau is approximately 0.18×10¹¹, 0.42×10¹¹ and 1.15×10¹¹ m³ yr⁻¹, respectively. From 1951 to 2017, the annual snowmelt increased significantly in the Tibetan Plateau and decreased significantly in North, Central and Southeast China. There was a decreasing trend in snowmelt in China, although this trend is not statistically significant. In West China, the mean annual snowmelt runoff ratio is more than 10 % in almost all third-level basins, except for those in the Taklimakan Desert. In basins in North and Northeast China, the snowmelt runoff ratio is generally more than 5 %, whereas in basins in South China it is less than 2 %. The Sen's slope shows that, from 1951 to 2017, the annual snowmelt runoff ratio decreased in most third-level basins in China; where this decrease is significant is mainly distributed in central Inner Mongolia, the southern slopes of the Tianshan mountains and South China.

Under RCP2.6, RCP4.5 and RCP8.5, the snowmelt in China shows a significant decreasing trend from 2006 to 2099, and the projected decrease in snowmelt in the near future (2011–2040; mid-future – 2041–2070; far future – 2071–2099) is 10.4 % (15.8 %; 13.9 %), 12.0 % (17.9 %; 21.1 %) and 11.7 % (24.8 %; 36.5 %), respectively, compared to the reference period (1981–2010). Among China's three main stable snow cover regions, the Tibetan Plateau is projected to experience the largest decrease in snowmelt in the future, followed by Northeast China; the snowmelt in the Northern Xinjiang is not projected to change significantly in the future. Under the three RCPs, the projected mean annual snowmelt runoff ratios in the third-level basins in the different future periods are less than those in the reference period, except for a few basins in Xinjiang and North China, with the largest decrease projected to occur in the Tibetan Plateau.

This study also investigated the spatial variability in snowmelt changes caused by changes in temperature and precipitation. It found that in low temperature regions (which can tolerate more warming), the snowmelt change was mainly influenced by precipitation, whereas in warm regions, the snowmelt change was most sensitive to temperature increases. The spatial variability in snowmelt changes may lead to the regional differences in the impact of snowmelt on water supply.

Although some uncertainties were introduced by the model principles, driving data and parameterization and other factors, this study is the first attempt to quantify snowmelt and its changes in China. Given the importance of snowmelt, the results have important implications for future research on snow–hydrology–climate interactions and contribute to water resource management planning under climate change.