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 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 10 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 pathways (RCP) scenarios (RCP2.6, RCP4.5 and RCP8.5). The results show that the mean annual snowmelt in China from 15 1951 to 2017 is 2.41×10 11 m 3 year -1 . The mean annual snowmelts in Northern Xinjiang, Northeast China, and the Tibetan Plateau – China’s three main stable snow cover regions – are 0.18×10 11 m 3 year -1 , 0.42×10 11 m 3 year -1 and 1.15×10 11 m 3 year -1 , respectively. From 1951 to 2017, the snowmelt increased significantly in the Tibetan Plateau and decreased significantly in North, Central and Southeast 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 20 in


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 phase change (Zhang and Ma, 2018;Pulliainen et al., 2020;Vano, 2020;You et al., 2020). Additionally, snow is a 35 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, 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 40 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 two types: physically-based 45 snowmelt models and simplified temperature index models (also known as degree-day methods) (Skaugen et al., 2018;. In theory, physically-based snowmelt models can provide more accurate predictions by considering the coupled interaction between energy components in complex snowmelt processes . 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., 50 2020). Temperature index models are based on the assumptions that the temporal variability of incoming solar radiation is adequately represented by the variations of 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., 55 2020).
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, Northeast China-Inner Mongolia (hereafter referred to as the 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 60 runoff in the Xinjiang Uygur Autonomous Region (hereafter, Xinjiang) is mainly supplied by snowmelt 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). Snow melting is also an important hydrological process on the Tibetan Plateau, which is the source region of many major Asian rivers and is considered as the Asian water towers (Immerzeel et al., 2010;Immerzeel et al., 2020). 65 Further, snowmelt is an important water source in other parts of China, especially the arid and semi-arid areas in North China 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;; however, country-wide assessments of snowmelt are lacking. Projected increases in air temperature and associated 70 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;Chen et al., 2019;Li et al., 2019;Zhang et al., 2020;Li et al., 2021), the spatiotemporal variability of 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 75 high spatial resolution (0.5′, approximately 1 km) dataset 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 as well as in its three main stable snow cover distribution regions; (2) analyse the spatiotemporal variability of snowmelt in China during 1951China during -2017 (3) elucidate the spatiotemporal variability of snowmelt runoff ratio in third-level basins in China; 80 (4) assess the impacts of projected future climate change on snowmelt in China.

Study region
In general, China can be divided into 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 85   (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 five types: prairie, taiga, tundra, mountain and ephemeral types   (Fig. 1c). The analysis of snowmelt runoff ratio is carried out at the three-level basin scale, 90 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).

High spatial resolution dataset of monthly air temperatures and precipitation 100
Data with a high spatial resolution of 0.5′ (approximately 1 km), including the monthly minimum, maximum, and mean temperatures (Tmin, Tmax, and Ta, respectively) and precipitation, were obtained from Zenodo in the Network Common Data across China, and the average values of the coefficient of determination (R 2 ), 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, respectively, and for the monthly precipitation are 0.863, 13.269 mm, 21.941 mm and 0.808, respectively. Detailed information on the dataset was given by Peng et al. (2019). Although the original data set covers the 110 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).

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/). 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. 120 The observational snow density data used to determine degree-day factors (DDF) were collected from 417 meteorological stations in China (Fig. 1c) during 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 125 interpolated using the parameters from the meteorological stations via inverse distance weighting (IDW) method.

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 130 (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  was obtained from the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn) (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 dataset 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 135 hydrological research in China (e.g. Liu et al., 2020;Wu et al., 2021;Zhu et al., 2021). The spatial resolution of this dataset is significantly different from that of the air temperature and precipitation data used in this study. The snow cover extent measures generated from this dataset was used to validate the model.
Additionally, a long time series of daily snow water equivalent dataset derived from passive microwave remote sensing data  was provided by the National Cryosphere Desert Data Center (https://www.ncdc.ac.cn) (Jiang et al., 2020) 140 ( Table 1). The snow water equivalent dataset was produced from the passive microwave remote sensing data using the mixed-pixel method. The dataset covers the entire land surface of China with a spatial resolution of 25 km. Detailed information on the dataset sources and product processes are given by Jiang et al. (2014) and Yang et al. (2019).

Climate projections and downscaling
Five CIMP5 models under three different representative concentration pathways scenarios (RCP2.6, RCP4.5 and RCP8.5), 145 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 (ISI-MIP) selected these five CMIP5 models to span the space of global mean temperature change and relative precipitation changes (Warszawski et al., 2014). Aforementioned climate projections have been bias-corrected downscaled to a grid with a resolution of 0.5° by ISI-MIP (Hempel et al., 2013), and widely utilized in climate change research in China (e.g. Li et al., 2016;Yuan et al., 2017;150 Chang et al., 2020). In this study, we use the delta downscaling method to determine the monthly future meteorological data  based on the high-spatial-resolution temperature and precipitation dataset and the simulations of the five CIMP5 models during the historical period . 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 its three main stable snow cover regions can be seen in Fig. S3 and Fig. S4, respectively. 155

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: 160 where M, Sacc, Psnow and Ssub are the snowmelt, snow accumulation, snowfall and snow sublimation (mm), respectively. DDF is the degree-day factor (mm ℃ -1 day -1 ), an empirical factor that relates the rate of snowmelt to air temperature. PDD is the accumulated positive air temperature (℃). The subscript m indicates the month, and D is the number of days in the month m.

Snowfall 165
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 mix, depending on two threshold temperature parameters: where Psnow and P are the monthly snowfall and precipitation (mm), respectively. Ta is the monthly mean temperature (℃). 170 Train and Tsnow are the threshold temperature (℃) for rainfall and snowfall, respectively.

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): where and are the density of snow and water (g cm -3 ) , respectively. The observed snow density is variable and the density of water is constant. DDF is the degree-day factor (mm ℃ -1 day -1 ). For taiga regions (Fig. 1c), the equation was corrected as follows: The DDF values at the meteorological stations were calculated for each month, and were then interpolated to each calculated 180 cell.

PDD
Based on the daily mean air temperatures from 1960 to 2018 at meteorological stations in China,  proposed a simple method to calculate the monthly PDD according to the monthly mean temperature: where T1 and T2 are threshold temperatures (℃); a, b and c are empirical parameters.
First, T1, T2, 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 190 measured monthly PDD, taking the MPZ as an example. 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 the Eq. 6 can adequately calculate the monthly PDD (Table 2).

Snow sublimation 205
As a loss of water from the snowpack to the atmosphere, snow sublimation is difficult to quantify by measurement or modelling (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. 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 210 sublimation to snowfall (e.g. Zhang et al., 2008;Zhu et al., 2014;Sexstone et al., 2018). Therefore, considering the limited data availability and monthly scale of this study, the following simple equation was used to estimate the monthly snow sublimation.
where k•Sacc and PET are the amount of snow accumulation available for sublimation and the potential sublimation (mm), 215 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.
PET was calculated by the Hargreaves-Samani method (Hargreaves and Samani, 1985): where Ta, Tmax and Tmin are the monthly average temperature, maximum temperature and minimum temperature (℃), 220 respectively; Ra, is the extraterrestrial radiation (MJ m -2 month -1 ), which is a function of the latitude and the Julian day; λ is the latent heat of vaporisation (MJ kg -1 ).

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 230 ratio, such as the ratios of total snowfall to total precipitation, total snowfall to total runoff, or melt season runoff to total annual runoff . 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: where M and Prain are the snowmelt and rainfall (mm) , respectively, and Snowr is the percentage of runoff derived from 235 snowmelt to the total runoff.

Trend analysis
The non-parametric Mann-Kendall test (Mann, 1945;Kendall, 1975) was used to analyse the trend and significance level of the snowmelt and other variables: where n is the number of data set, xj and xk are the data values in time series j and k, m is the number of tied groups and ti denotes the number of ties of extent i. A tied group is a set of sample data having the same value. Positive and negative 245 values of Z indicate increasing and decreasing trends, respectively. Testing trends is done at the specific α significance level.
If | | > 1− /2 , 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 slope method (Sen, 1968) was used to analyse the slope of the variation: where β sign reflects whether a trend is negative or positive, while its value indicates the steepness of the trend.

Evaluation criteria
Statistical indices were used for quantitative analysis of the snowmelt model performance, and a series of statistical criteria used in this study as follows: where R 2 , MAE, RMSE and NSE are the coefficient of determination, mean absolute error, root mean square error, and Nash-Sutcliffe efficiency, respectively. Xsi and Xoi represent simulated and observed data at time i, respectively. n is the number of 260 data points, and cov and σ are the covariance and standard deviation, respectively.

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. Fig. 4 shows the statistical analyses of the calculated and observed snow depths in China from 1951 to 2009. Similar to snowfall validation, besides R 2 and NSE, the ratio of MAE and RMSE to the mean snow depth 280 (MAE/Mean and RMSE/Mean, respectively) were selected to analyse the model performance.

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

Snow water equivalent validation using the dataset 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 dataset 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 m 3 for different regions. The mean monthly snow accumulation in China output by the model from 1980 to 2017 is 305 2.55×10 10 m 3 , while the mean snow water equivalent derived from passive microwave remote on the last day of each month is 1.47×10 10 m 3 . The R 2 , MAE, RMSE and NSE between the snow accumulation output by the model against snow water equivalent dataset in China from 1980 to 2017, are 0.62, 1.27×10 10 m 3 , 1.27×10 10 m 3 and 0.80, respectively. The R 2 , 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 dataset 310 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 modelled 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 315 whole of China (MAE and RMSE decreased and NSE increased; see Fig. S7). respectively. Southeast China (mainly the Huaihe River basin and the lower reaches of the Yangtze River) also has high 325 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.
The spatial pattern of snowmelt differs considerably during the year (Fig. 8)  The month of maximum snowmelt differs between regions (Fig. 7b). In the 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 340 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 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 345 elevations mountains in the Tibetan Plateau, with air temperatures above 0 ℃ for only a few days in the summer, and the month of maximum snowmelt in those high elevations areas also occurs in summer.

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 355 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 (Fig. 7c, Fig. 7d and Fig. 9). From 1951 to 2017, the linear trends of annual snowmelt in China, Northern Xinjiang, Northeast China, and Tibetan Plateau are -2.7×10 9 m 3 decade -1 , 0.2×10 9 m 3 decade -1 , -0.3×10 9 m 3 decade -1 and 1.4×10 9 m 3 decade -1 , respectively. 360

Figure 9. Interannual variability of the mean annual snowmelt in China and its three main stable snow cover regions from 1951 to 2017.
The temporal trend of snowmelt in each month from 1951 to 2017 varies across China (Fig. S9, Fig. S10). In Northern Xinjiang and Northeast China, snowmelt increased significantly in March, but decreased significantly in April. This might 365 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.

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 region -Northern Xinjiang and the Tibetan Plateauthe mean annual snowmelt runoff ratios are more than 10% in all basins, except those in the Taklimakan Desert. In the basins in North and 375 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).
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 cold months, and there is no 380 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 385 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.

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 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 400 (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 Gurbantunggut 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 4 basins had significant increasing trends: the source region of the Yellow River and three sub-basins of the Songhua River.
The temporal trend in the snowmelt runoff ratio from 1951 to 2017 shows spatial variations in China in every month (Fig.  405   S12, Fig. 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 sub-basins 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 the 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 410 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 Mountain and 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 415 runoff ratios showed a large increase in the upper reaches of four rivers: 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).  Table 4. In Northern 425

Future changes of snowmelt
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 . 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 snowmelts in different future periods under the three RCPs are all 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 430 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%), 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 435 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.

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 455 general, the largest decreases in snowmelt runoff ratio in the basins are projected to occur by the 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 460 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).

Model Evaluation 470
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 2 values between the simulated snowfall and the snowfall observed at three 475 meteorological stations of 0.42, 0.34 and 0.61, respectively; in this study, the R 2 values at the same stations are 0.72, 0.46 and 0.39, respectively. Li et al. (2020) used temperature thresholds to calculate the snowfall in the Tianshan Mountains of Central Asia, and obtained a mean R 2 value between the simulated snowfall and the snowfall observed at 27 meteorological stations of 0.61; in this study, the mean R 2 value at 50 stations in Xingjiang 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 2 values 480 between simulated and observed snowfall of less than 0.3 for most meteorological stations; in this study, the R 2 is larger than 0.3 for most stations. Snow depth is the worst-performing model output, 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 485 scale, where the snow properties are not always representative of snow at grid (Sexstone et al., 2020). (3) Snow depth itself is difficult to simulate, and even the snow depth retrieved by remote sensing have 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 cm 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 AMSR-E passive 490 microwave against ground observations from meteorological stations across China, and found that the RMSE varied from 6.85 cm to 16.79 cm in Xinjiang region, and from 6.21 cm to 18.05 cm in Northeast China. In this study, the RMSE varies from 0.56 cm to 13.32 cm in Xinjiang, and 0.54 cm 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 495 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 2 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 500 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 505 is reliable for calculating the snowmelt in China.

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 510 2017 were selected to analyse 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 girds where snowfall increased significantly and snowmelt decreased significantly. In a few grid cells, snowfall decreased significantly and the 515 snowmelt increased significantly. These grids are too few and can be ignored (Table S2). Compared with the total annual precipitation, the change of snowfall has more influence on snowmelt. Tan et al. (2019) analysed 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 520 annual scale (e.g. Woldemeskel et al., 2016;Ahmed et al., 2019;Yue et al., 2019). Assessing precipitation from climate models using multiple time scales or seasons may allow a more accurate prediction of snowfall and snowmelt changes.  Table S2.

525
As snow cover formation and snowmelt are closely related to the temperature threshold of 0 ℃, temperature change near 0 ℃ is more likely to trigger a drastic change in snowmelt. This may partly explain for 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 the freezing point (Qin et al., 530 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 535 under RCP2.6, RCP4.5 and RCP8.5. The climate in this region is warm with average monthly temperatures near 0 ℃ 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 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 540 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 CIMP5 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 545 in less snowfall and snowmelt at medium and low elevations, and the temperature may become a more influential control on snowmelt.

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 550 (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.

Northern Xinjiang
Northern Xinjiang is located in the TCZ, and snowmelt is an important source of fresh water in this region 555 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%, also show a significant decreasing trend. March and April are 560 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.

Northeast China
In Northeast China, the runoff in March and April is mainly generated from snowmelt. From 1951From to 2017 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 snow melting releases water in spring, which can influence agriculture by affecting soil moisture (Qi et al., 2020). Snowmelt is vital for seasonal water supply in Northeast China , which is an important region for agricultural production in China. 570 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.

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 575 Tibetan Plateau are the source of many arid inland rivers such as the Hetian 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 580 runoff. The snowmelt in the Tibetan Plateau showed an increasing trend from 1951 to 2017, and as the temperature is projected to rises 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 585 supplies.

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 590 caused a significant decrease in snowmelt in this region, and this trend is projected to continue in the future, the impact on water supply is likely to be much smaller than in other regions of China. The increased frequency and intensity of flooding that are projected to be caused by climate change is a major concern for this region (Qin and Lu, 2014).

Uncertainties and limitations
The snowmelt simulation reported in this study was largely driven by the climatology dataset of WorldClim using delta 595 spatial downscaling, and any uncertainties in the dataset 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 dataset generated by bilinear interpolation method containing detailed topographic information, as well as the effects of distance to the nearest coast and satellite-derived covariates, is of high quality, there is still a gap 600 between the downscaled data and the observed data, especially in the mountainous areas (Peng et al., 2019;Ding and Peng, 2020). Uncertainties in the CIMP5 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: the threshold temperatures for 605 the separation of rainfall and snowfall (Train and Tsnow), the degree-day factors (DDF) and the ratio of snow sublimation to snow accumulation (k). Although Train and Tsnow 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 610 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 modelling (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 sublimation.
However, this ratio can vary considerably both spatially and temporally, having been estimated to vary between 0.1 and 90% 615 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. 620 In addition to snowmelt, the snowmelt runoff ratio is another important research object in this study. We estimated this ratio using snowmelt / (snowmelt + rainfall). However, this method may underestimate the ratio because snowmelt is more effective at generating catchment runoff than rainfall Jenicek and Ledvinka, 2020).

Conclusions
A simple temperature index model was used to calculate a time series of monthly snowmelt at a resolution of 0.5′ in China 625 for the 1951-2017 period. The mean annual snowmelt in China is 2.41×10 11 m 3 year -1 , and the mean annual snowmelt in Northern Xinjiang, Northeast China and Tibetan Plateau is approximately 0.18×10 11 m 3 year -1 , 0.42×10 11 m 3 year -1 and 1.15×10 11 m 3 year -1 , respectively. From 1951 to 2017, the annual snowmelt increased significantly in the Tibetan Plateau, and decreased significantly in the 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 630 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 (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), respectively. 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 640 projected to change significantly in the future. Under the three RCPs, the projected mean annual snowmelt runoff ratios in the third-level basins in 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 of 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 645 precipitation, whereas in warm regions, snowmelt change was most sensitive to temperature increases. The spatial variability of 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 researches on snow-hydrology-climate interactions and contribute to water resource 650 management planning under climate change.
Data availability. The spatial resolution with 0.5′, including monthly minimum, maximum, and mean temperatures (Tmin, Tmax, and Ta) and precipitation, are obtained from Zenodo in Network Common Data Form at https://doi.org/10.5281/zenodo.3114194 for precipitation and https://doi.org/10.5281/zenodo.3185722 for air temperatures. 655 The air temperature, snowfall, snow depth and snow density observation data are obtained from the China Meteorological Administration (http://data.cma.cn/). The snow depth dataset derived from passive microwave remote sensing data is obtained by the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn). The snow water equivalent dataset derived from passive microwave remote sensing data is obtained by the National Cryosphere Desert Data Center (https://www.ncdc.ac.cn). The CMIP5 data set is distributed by the Inter-Sectoral Impact Model Intercomparison Project on 660 its own website (http://www.isi-mip.org). The third-level basins dataset is obtained by the Resource and Environment Science and Data Center (https://www.resdc.cn/).
Author contributions. RC initiated the study. YY, RC and GL developed the methodology. YY performed all analyses. YY and ZL prepared the input meteorological data. GL and XW prepared the CIMP5 models data. YY prepared the manuscript with contributions from RC. 665 Competing interests. The authors declare that they have no conflict of interest.